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Abstract 

This  thesis  used  vortex  lattiee  lifting  line  theory  to  model  an  axisymmetrical-dueted  propeller 
with  no  gap  between  the  duet  and  the  propeller.  The  theory  required  to  model  the  duet  and  its 
interaetion  with  the  propeller  were  discussed  and  implemented  in  Open-source  Propeller  Design 
and  Analysis  Program  (OpenProp).  Two  routines  for  determining  the  optimum  circulation 
distribution  were  considered,  and  a  method  based  on  calculus  of  variations  was  selected.  The 
results  of  this  model  were  compared  with  the  MIT  Propeller  Eifting  Fine  Program  (PEE)  output 
for  the  purpose  of  validation. 

Ducted  propellers  are  prevalent  in  modem  marine  propulsion  systems,  and  the  application  of  this 
technology  continues  to  expand.  The  theory  associated  with  ducted  propellers  applies  to  a  wide- 
range  of  devices  which  include  azimuth  thmsters,  pumpjets,  and  tidal  turbines.  Regardless  of  the 
application,  engineers  need  tools  such  as  OpenProp  to  design  these  devices  for  their  expected 
operating  conditions.  OpenProp  is  an  open  source  MATLAB®-based  suite  of  propeller 
numerical  design  tools.  Previously,  the  program  only  designed  open  propellers.  The  code 
developed  in  this  thesis  extended  OpenProp’s  capability  to  be  able  to  design  a  propeller  within 
an  axisymmetrical  duct. 


Thesis  Supervisor;  Patrick  J.  Keenan 
Title:  Professor  of  Naval  Architecture 

Thesis  Supervisor:  Richard  W.  Kimball 


2 


Acknowledgements 


The  author  thanks  the  following  individuals  for  all  of  their  support  and  assistance  with  this 
thesis: 

Professor  Rich  Kimball  for  all  his  wisdom  and  guidance  during  not  only  the  thesis  process  but 
also  during  the  two  classes  he  taught.  His  classes  and  teaching  method  were  among  the  best 
experienced  at  MIT. 

CAPT  Patrick  Keenan  for  the  leadership  and  direction  he  provided  both  during  the  thesis  process 
and  throughout  the  entire  course  2N  program.  His  course  was  also  one  of  the  best  the  author 
experienced  at  MIT. 

Brenden  Epps  for  his  assistance  with  the  MATLAB®  coding  involved  with  modeling  ducted 
propellers. 

And  most  of  all,  my  family  for  all  of  their  support  and  understanding  during  my  time  at  MIT. 
While  Christine  and  Jake  experienced  the  entire  adventure,  Hank  arrived  just  in  time  to  see  the 
thesis  completed. 


3 


3 


Table  of  Contents 

Acknowledgements .... 

List  of  Figures . 6 

1 .  Introduction . 8 

2.  Overview  of  the  Propeller  Design  code:  OpenProp . 1 1 

3.  Theoretical  Foundation . 16 

3.1  Ducted  Propeller  Theory . 16 

3.2  Vortex  Ring  Theory  and  Algorithm . 22 

3.3  Circumferential  Mean  Veloeity . 25 

3.4  Cireulation  Optimization . 32 

4.  Implementation  and  Validation . 37 

5.  Conelusions  and  Reeommendations . 47 

5.1  Conclusions . 47 

5.2  Recommendations  for  further  work . 48 

References . 49 

Appendix  A.  Duct  Theory  MATLAB®  Code . 51 

A.l  ductVort.m . 51 

A.2  ductThrust.m . 53 

A. 3  vRing.m . 55 

A.4  vpfDuct.m . 57 

A.5  CMV.m . 58 

A.  6  ductPlot.m . 60 

Appendix  B.  Mathematical  Functions  MATLAB®  Code . 63 

B. l  Q2half.m . 63 

B.2  Q2Mhalfm . 63 

4 


B.3  Heuman.m 


64 


Appendix  C.  Variational  Optimization  Routine  MATLAB®  Code . 65 

C.l  Coney. m . 65 

C.2  Align_wake.m . 71 

Appendix  D.  Test  Case  Setup  Data . 72 


5 


List  of  Figures 

Figure  1-1:  Ducted  propeller  and  associated  vortex  system  representation  from  Coney  (1) . 9 

Figure  2-1:  OpenProp’s  parametric  analysis  input  GUI . 12 

Figure  2-2:  Efficiency  diagrams  produced  by  OpenProp’s  parametric  analysis . 12 

Figure  2-3:  OpenProp’s  single  propeller  design  input  GUI . 13 

Figure  2-4:  OpenProp’s  graphical  out  of  key  propeller  parameters . 14 

Figure  2-5:  Blade  and  propeller  representations  from  OpenProp . 14 

Figure  3-1 :  Diagram  of  the  ducted  propeller  vortex  system  from  Coney  (1) . 17 

Figure  3-2:  Optimum  circulation  distributions  for  ducted  and  open  propellers . 19 

Figure  3-3:  Optimization  circulation  distribution  for  a  zero  gap  ducted  propeller . 20 

Figure  3-4:  CMV.m  Validation  Representative  Circulation . 28 

Figure  3-5:  Heuman’s  Lambda  Function . 31 

Figure  3-6:  OpenProp  t=0.80  test  case  using  the  Lerbs-based  optimization  routine . 33 

Figure  4-1:  OpenProp  v2  single  propeller  design  GUI  with  duct  parameters . 37 

Figure  4-2:  OpenProp  v2  required  duct  parameters . 37 

Figure  4-3:  Sample  rendering  of  ducted  propeller  produced  in  the  test  cases . 39 

Figure  4-4:  OpenProp  v2  algorithm  propeller  results  using  PEL  circulation . 40 

Eigure  4-5:  OpenProp  v2  algorithm  duct  results  using  PEL  circulation . 41 

Eigure  4-6:  Efficiency  versus  thrust  ratio  (x)  comparison  between  OpenProp  v2  and  PEL . 42 

Eigure  4-7:  OpenProp  v2  graphical  output  for  test  case  with  t  =  0.8 . 43 

Eigure  4-8:  OpenProp  v2  comparison  with  PEL  for  t  =  0.8 . 43 

Eigure  4-9:  OpenProp  v2  comparison  with  PEL  for  t  =  0.8  (duct  ring  velocities) . 44 

Eigure  4-10:  OpenProp  v2  comparison  with  PEL  for  t  =  1.0 . 45 

Eigure  4-11:  OpenProp  v2  comparison  with  PEL  for  t  =  1.2 . 46 

Eigure  D-1:  Test  case  parameters . 72 

Eigure  D-2:  OpenProp  v2  input  GUI  for  test  cases . 72 

Eigure  D-3:  PEL  current  settings  for  inviscid  and  viscid  test  cases . 73 


6 


Figure  D-4:  PLL  overall  input  file  for  test  eases . 73 

Figure  D-5:  Sample  PLL  output  summary  for  test  case  run . 74 


7 


1.  Introduction 


Ducted  propellers  are  widely  used  in  marine  propulsion  systems  for  a  variety  of  reasons.  As 
shown  by  Kort  in  1934,  dueled  propellers  ean  aehieve  higher  efficieney  partieularly  in  slow, 
heavily  loaded  applieations  sueh  as  tugboats  and  oeean  platforms.  While  this  effieieney  gain  is 
lost  at  higher  speeds  due  to  the  viseous  drag  from  the  duet,  there  are  many  more  reasons  that  a 
designer  might  ehoose  a  dueled  propeller  or  derivative  sueh  as  a  pump  jet  or  water  jet  instead  of 
a  traditional  open  propeller.  Several  of  these  reasons  are  listed  below. 

Greater  power  density.  A  dueled  propeller  ean  produee  more  thrust  than  an  open 
propeller  of  the  same  size.  If  ship  geometry  limits  the  size  of  the  propulsor,  a  dueled 
propeller  might  be  the  best  option. 

Physieal  proteetion.  A  duet  provides  proteetion  for  the  propeller  blades. 

Cavitation  reduetion.  As  deseribed  below,  a  deeelerating  duet  ean  be  used  to  inerease 
the  statie  pressure  at  the  propeller  whieh  reduees  or  eliminates  eavitation. 

Simplieity.  Beeause  dueled  propulsors  often  ineorporate  direetional  eontrol  (veetored 
thrust  via  a  trainable  duet,  steerable  nozzle,  eet),  rudders  ean  be  eliminated. 

Expanded  operational  environment.  Beeause  the  duet  provides  proteetion  and  ean 
reduee  or  eliminate  other  appendages  sueh  as  a  rudder,  dueled  propellers  ean  improve 
shallow  water  operation.  This  is  espeeially  true  for  a  water  jet  sinee  only  the  inlet 
must  be  submerged  for  proper  operation. 

Another  promising  use  of  dueled  propeller  teehnology  is  the  tidal  turbine  market.  With  the 
world  looking  inereasingly  towards  renewable  energy  sourees,  harnessing  the  power  of  the  oeean 
is  of  great  interest  and  importanee.  Development  of  tidal  turbines  direetly  leverages  the  researeh 
assoeiated  with  dueled  propellers.  The  geometries  of  both  problems  are  essentially  the  same,  and 
the  design  proeess  of  eaeh  involve  an  optimization  involving  thrust  and  torque.  For  dueled 
propulsors,  the  goal  is  provide  a  eertain  amount  of  thrust  while  minimizing  the  required  torque  a 
ship’s  engines  must  provide  through  a  shaft  or  eleetrie  motor.  The  optimization  is  essentially 
reversed  for  a  tidal  turbine  that  extraets  power  from  the  oeean  via  a  turbine  generator.  In  this 
ease,  the  designer  desires  to  maximize  torque  and  minimize  thrust.  While  this  thesis  develops 
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the  design  tool  for  a  traditional  dueled  propeller,  the  eoneepts  and  most  of  the  eode  are  direetly 
applieable  to  tidal  turbines. 

Ducted  propellers  and  associated  derivatives  (electric  drive,  pods,  azimuth  thrusters,  water  jets, 
pumpjets,  etc)  will  continue  to  play  an  important  role  in  ship  propulsion  and  will  have  an 
expanding  role  in  renewable  energy  efforts  as  described  above.  For  this  reason,  engineers  need 
tools  to  design  the  devices  for  specific  operating  conditions. 

Following  the  approach  presented  by  Coney  in  (1),  this  thesis  developed  the  MATLAB® 
algorithms  necessary  to  model  a  duct  and  its  interactions  with  a  propeller.  These  algorithms 
were  then  integrated  into  the  existing  OpenProp  propeller  design  program.  The  model  assumed 
that  there  was  no  gap  between  the  duct  and  propeller.  The  duct  was  represented  by  an  image 
system  of  vorticity  and  a  system  of  ring  vorticies  at  the  radius  of  the  duct  cylinder  (Figure  1-1). 


-c  -  c. 


-r(M)c> 


Image  Vortex  System 


Ring  Vortex  System 


£7 


Figure  1-1:  Ducted  propeller  aud  associated  vortex  system  represeutatiou  from  Couey  (1) 
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The  image  system  modeled  the  nonaxisymmetric  effeet  of  the  duct  while  the  ring  vorticies 
provided  an  estimate  of  the  resulting  duct  force.  The  influence  functions  calculated  for  the  radial 
lifting  line  control  points  included  the  effects  from  the  duct  image  system,  and  the  inflow  was 
modified  by  the  effect  of  the  duct  ring  vorticies.  A  variational  optimization  routine  was 
employed  to  determine  the  optimum  circulation  distribution  for  the  lifting  line.  The  model 
accounted  for  viscous  drag  but  duct  thickness  was  neglected.  To  the  greatest  extent  possible, 
this  thesis  used  the  notation  presented  in  (1). 
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2.  Overview  of  the  Propeller  Design  code:  OpenProp 


Open-source  Propeller  Design  and  Analysis  Program  (OpenProp)^  is  an  open  source 
MATLAB®-based  suite  of  propeller  numerical  design  tools.  This  program  is  an  enhanced 
version  of  the  MIT  Propeller  Vortex  Lattice  Lifting  Line  Program  (PVL)  developed  by  Professor 
Justin  Kerwin  at  MIT  in  2001 .  OpenProp  vl  .0,  originally  titled  MPVL,  was  written  in  2007  by 
Hsin-Lung  Chung  and  Kate  D’Epagnier  and  is  described  in  detail  in  (2)  and  (3).  Two  of  its  main 
improvements  versus  PVL  are  its  intuitive  graphical  user  interfaces  (GUIs)  and  greatly  improved 
data  visualization  which  includes  graphic  output  and  three-dimensional  renderings. 

OpenProp  was  designed  to  perform  two  primary  tasks:  parametric  analysis  and  single  propeller 
design.  Both  tasks  begin  with  a  desired  operating  condition  defined  primarily  by  the  required 
thrust,  ship  speed,  and  inflow  profile.  The  parametric  analysis  produces  efficiency  diagrams  for 
all  possible  combinations  of  number  of  blades,  propeller  speed,  and  propeller  diameter  for  ranges 
and  increments  entered  by  the  user.  The  efficiency  diagrams  are  then  used  to  determine  the 
optimum  propeller  parameters  for  the  desired  operating  conditions  given  any  constraints  (e.g. 
propeller  speed  or  diameter)  specified  by  the  user.  Figure  2-1  shows  the  input  GUI  for  the 
parametric  analysis  routine,  and  Figure  2-2  shows  the  efficiency  diagrams  produced  by  that 
routine. 


*  Throughout  this  thesis,  OpenProp  refers  to  the  design  program  in  general.  OpenProp  vl.O  (version  1.0)  refers  to 
the  original  version  of  OpenProp  which  was  developed  for  open  propellers  only,  and  OpenProp  v2.0  refers  to  the 
version  associated  with  this  thesis  which  includes  the  capability  to  model  ducted  propellers. 
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Figure  2-1:  OpenProp’s  parametric  analysis  input  GUI 


Number  of  Blades:  3 


Number  of  Blades:  4 


Number  of  Blades:  6 


Figure  2-2:  Efficiency  diagrams  produced  by  OpenProp’s  parametric  analysis 
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The  single  propeller  design  routine  produces  a  complete  propeller  design  for  the  desired 
operating  condition  and  defined  propeller  parameters  (number  of  blades,  propeller  speed, 
propeller  diameter,  hub  diameter,  etc).  Figure  2-3  shows  the  input  GUI  for  the  single  propeller 
design  routine.  OpenProp’s  graphical  out  of  key  propeller  parametersFigure  2-4  shows  the 
graphical  output  of  key  propeller  parameters,  and  Figure  2-5  shows  blade  profiles  and  complete 
three-dimensional  representation  of  the  propeller. 
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Figure  2-3:  OpenProp’s  single  propeller  design  input  GUI 
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J=0.750;  Ct=0.994;  Cq=0.346;  Kt=0.220;  Kq=0.038:  ti  =0.686 
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Figure  2-4:  OpenProp’s  graphical  out  of  key  propeller  parameters 
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OpenProp  was  developed  to  serve  as  an  open  souree  eode  for  propeller  design.  While  it  is 
currently  a  tool  that  is  only  used  in  the  initial  design  phase,  it  is  a  base  program  that  can  be 
continually  expanded  to  perform  detailed  design  and  analysis  of  sophisticated  marine  propulsors 
and  turbines.  Extending  OpenProp  to  include  a  duct  was  the  main  focus  of  this  thesis. 
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3.  Theoretical  Foundation 


3.1  Ducted  Propeller  Theory 

Ducted  propellers  are  generally  divided  into  two  types:  aeeelerating  and  deeelerating  nozzles. 
The  aeeelerating  duet  or  “Kort”  nozzle  has  been  widely  used  sinee  Kort  showed  in  1934  that  this 
type  of  duet  produees  a  positive  thrust  and  ean  inerease  effieieney  in  heavily  loaded  applieations 
sueh  as  tugboats.  It  was  also  shown  that  the  optimum  diameter  for  a  dueted  propeller  is  smaller 
than  that  for  an  open  propeller.  Beeause  of  this,  aeeelerating  duets  are  sometimes  used  when 
inereased  thrust  is  needed  from  a  propeller  whose  size  is  eonstrained  by  the  ship’s  eharaeteristies 
or  operating  eonditions. 

A  deeelerating  duet  inereases  the  statie  pressure  at  the  propeller  and  is  used  to  reduee  eavitation. 
A  reduetion  in  eavitation  lowers  the  noise  generated  by  a  propeller  and  reduees  erosion  of  the 
blades. 


From  momentum  theory,  the  ideal  effieieney,  iji,  of  a  dueted  propeller  is  given  in  Equation  3-1 
where  t  is  the  thrust  ratio  (Equation  3-2),  Cj  is  the  thrust  eoeffieient  (Equation  3-3),  p  is  the  fluid 
density.  Vs  is  the  ship  speed,  and  D  is  the  propeller  diameter  (4). 


2 

1  -f  1  -f  tCj 

3-1 

Tp  Propeller  Thrust 
T  Total  Thrust 

3-2 

T 

-  T  - 

3-3 
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As  the  thrust  ratio  is  lowered,  the  duct  produces  more  of  the  required  total  thrust  and  the  ideal 
efficiency  increases.  However,  when  finite  blade  effects  are  considered  a  penalty  is  paid  for  the 
increased  axial  velocity  at  the  propeller  plane  and  efficiency  decreases  after  reaching  a  maximum 
at  approximately  t  =  0.9  (1).  Additionally,  since  the  thrust  ratio  is  multiplied  by  the  thrust 
coefficient,  a  large  thrust  coefficient  is  required  in  order  to  realize  a  significant  efficiency  gain. 
Hence,  ducted  propellers  are  commonly  used  in  heavily  loaded  situations  such  as  tugboats. 

As  introduced  above,  this  thesis  modeled  the  duct  as  an  infinite  cylinder  by  adding  an  image 
vortex  system  to  the  vortex  lattice  representing  the  lifting  line  and  adding  a  system  of  ring 
vortices  to  account  for  duct  forces  and  the  axisymmetric  mean  inflow  modification  by  the  duct. 
Figure  3-1  shows  the  ducted  propeller  vortex  system  used  by  Coney  and  implemented  in 
OpenProp.  The  propeller  coordinate  system  is  shown  in  the  lower  left-hand  corner  of  Figure 
3-1. 


Figure  3-1:  Diagram  of  the  ducted  propeller  vortex  system  from  Couey  (1) 


^  The  propeller  coordinate  system  used  in  OpenProp  and  this  thesis  has  positive  x  in  the  downstream  direction. 
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The  image  system  method  was  essentially  the  same  as  the  method  used  originally  in  OpenProp  to 
model  the  hub.  The  underlying  assumption  here  was  that  the  image  system  will  approximately 
satisfy  the  eondition  of  zero  radial  veloeity  at  the  eylinder.  The  image  system  used  a  constant 
pitch  angle  based  on  the  tip  trailing  vortex  on  the  lifting  line. 

The  radii  of  the  image  vortices  and  their  associated  tan()ffj)  were  determined  using  Equations  3-4 
and  3-5  where  is  the  radius  of  the  duct  image  shed  vortex  trailer,  is  the  radius  of  the  duct 

cylinder,  r  is  the  radius  of  the  helical  trailing  vortex  shed  by  the  lifting  line,  Pi  is  the 
hydrodynamic  pitch  angle,  and  the  subscript  ‘"tip"  refers  to  the  or  last  shed  vortex  trailer  on 
the  lifting  line. 


3-4 


=  tan  (Pi) tip 


3-5 


The  radial  and  tangential  influence  functions  from  the  duct  image  were  added  to  the  radial  lifting 
line  influence  functions  as  follows: 

total  Xl(i(ji,Tlx)  -f  duct  image 

3-6 


\Ut(TL,Tyi)^total  Ut(Tl,TTl)  -f  duct  image 


3-7 


These  influence  functions  were  used  in  the  variational  optimization  routine  described  below  to 
determine  the  optimum  circulation  distribution  for  ducted  propellers. 

As  the  gap  between  the  duct  and  propeller  tip  is  decreased,  the  optimum  circulation  distribution 
becomes  more  tip  loaded.  At  the  limiting  case  of  zero  tip  gap,  the  circulation  reaches  its 
maximum  value  at  the  tip.  Figure  3-2  shows  the  OpenProp  v2  and  PEL  results  for  the  optimum 


18 


circulation  distribution  for  neutrally  loaded  (t  =  1.0)  ducted  propellers  (Cj  =  1.2,  Js  —  0.6) 
represented  by  a  duct  image  system.  The  OpenProp  v2  results  were  essentially  identical  to  those 
obtained  from  PLL. 


Figure  3-2:  Optimum  circulatiou  distributious  for  ducted  and  open  propellers 

In  (1),  Coney  compared  the  duct  image  method  with  a  more  sophisticated  panel  method 
representation  of  the  duct  and  determined  that  the  optimum  circulation  distribution  agreed  very 
well.  Since  a  panel  method  would  be  too  computationally  intensive  for  an  early  stage  design  tool 
such  as  OpenProp,  the  results  were  extremely  fortuitous.  Coney’s  results  (1)  are  shown  in  Figure 
3-3. 


19 


m 
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Figure  3-3:  Optimization  circulation  distribution  for  a  zero  gap  ducted  propeller  as  determined  from  image 

and  panel  representations  of  the  duct  from  Coney  (1). 

The  ring  vortex  system  was  used  to  estimate  the  duet  foree  and  the  mean  modifieation  of  the 
inflow  at  the  propeller  lifting  line.  The  length  of  the  system  represented  the  duct  chord  length, 

Crf,  which  was  chosen  to  be  equal  to  the  propeller  radius,  R.  The  ring  vortices  of  the  system  were 
spaced  (As)  evenly  at  approximately  the  same  constant  interval  used  to  discretize  the  lifting  line. 
Furthermore,  the  system  was  positioned  such  that  the  lifting  line  was  located  at  the  duct  mid¬ 
chord  and  between  ring  vortices.  The  strengths  of  ring  vorticies  were  calculated  to  represent  a 
NACA  a=0.8  meanline  over  the  length  of  the  duct. 


The  axial  component  of  the  velocity  that  the  ring  vortex  system  induced  on  the  lifting  line  is 
shown  in  Equation  3-8  where  is  the  number  of  ring  vortices  used  to  represent  the  duct.  Fd  is 
the  strength  of  the  ring  vortex,  n)  is  the  velocity  induced  by  the  vortex  ring  of 

unit  strength  on  the  propeller  lifting  line  control  point.  The  algorithm  used  to  calculate  the 
induced  velocities  from  the  vortex  rings  is  discussed  in  the  next  section. 
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3-8 


Nd 

^apd^O  =  ^  Kpd  ii-n)Td,(n) 

n=l 


The  axial  force  on  each  vortex  ring  was  calculating  by  using  the  Kutta-Joukowski  law.  As 
shown  in  Equation  3-9,  the  total  thrust,  on  the  duct  was  calculated  by  summing  the  axial 
forces  on  the  vortex  rings  and  adding  a  viscous  drag  by  using  a  two-dimensional  airfoil  sectional 
drag  coefficient, 


Nd 

Ta  =  -2pnra  +  Urin,raWd(n)  +  [Va(ra)  +  Ua(n,ra)]Cd,CQj  As 

n=l 

3-9 


Here,  14- (?d)  ^^e  the  radial  and  axial  inflow  velocities,  respectively,  at  the  duct 

radius,  r^.  Uy(n,  and  Ua(n,  are  the  radial  and  axial  components  of  the  circumferential 
mean  velocity  induced  on  the  ring  vortex  by  the  propeller.  Calculation  of  the  circumferential 
mean  velocities  is  discussed  below. 

The  theory  above  was  implemented  in  OpenProp  v2  using  two  new  MATLAB®  functions, 
ductVort.m  and  ductThrust.m  (Appendix  A).  ductVort.m  was  used  to  determine  the  vortex  ring 
system  that  represented  the  duct  and  to  calculate  the  induced  velocity  on  the  lifting  line  from  that 
system.  ductThrust.m  was  used  to  calculate  the  total  duct  thrust  coefficient  derived  from 
equation  3-9  that  was  required  in  the  optimization  routine  discussed  below.  It  also  scaled  the 
duct  vortex  ring  system  circulation  so  that  the  duct  provides  the  thrust  specified  by  the  thrust 
ratio. 
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3.2  Vortex  Ring  Theory  and  Algorithm 


OpenProp  uses  vortex  rings  to  model  the  duet’s  axisymmetric  mean  modifieation  of  the  inflow 
to  the  propeller  and  to  ealeulate  the  thrust  provided  by  the  duet.  A  vortex  ring  is  a  vortex 
filament  that  forms  a  eirele  of  radius  r'.  Applying  the  Biot-Savart  law,  Kuehemann  and  Weber 
(5)  derived  the  infiuenee  from  a  vortex  ring  in  terms  of  eomplete  elliptie  integrals,  K(k)  and  E(k). 
K(k)  denotes  the  eomplete  elliptie  integral  of  the  first  kind,  and  E(k)  denotes  the  eomplete  elliptie 
integral  of  the  seeond  kind.  The  argument,  k,  is  known  as  the  modulus. 

As  derived,  the  vortex  ring  is  in  the  Y-Z  plane  and  is  loeated  at  x  =  0.  Due  to  the  inherent 
symmetry  of  a  ring,  all  points  contained  on  a  circle  in  the  Y-Z  have  the  same  induced  velocity. 
Therefore,  Kuehemann  and  Weber  derived  equations  for  the  axial  and  radial  components  of  a 
vortex  ring  with  radius  r'and  strength  F  for  a  point  located  on  a  circle  of  radius  i?  in  a  plane  X 
units  from  and  parallel  to  the  Y-Z  plane.  Additionally,  the  center  of  the  circle  is  on  the  x-axis. 
The  axial  induced  velocity  is  given  in  Equation  3-10  with  the  arguments  and  variables  shown  in 
Equations  3-11  and  3-12.  The  radial  induced  velocity  is  given  in  Equation  3-13. 


Vy^ix,  r) 


r  1 
27rr' ^^2  +  (r  +  1)2 


1  + 


2(r  -  1) 
x2  -I-  (r  —  1)2 
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For  OpenProp,  these  equations  were  implemented  in  a  function  named  vRing.m  (Appendix  A). 
vRing.m  was  successfully  validated  using  Tables  1  and  2  contained  in  (6).  vpfDuct.m  (Appendix 
A)  is  a  velocity  prediction  function  for  a  duct  modeled  only  with  vortex  rings  (i.e.  no  thickness). 
A  future  improvement  will  include  source  rings  to  model  duct  thickness. 

MATLAB®  contains  a  defined  function,  ELLIPKE,  for  calculating  the  elliptic  integrals  of  the 
first  and  second  kind.  However,  the  argument  for  this  function  is  not  the  modulus,  k.  Instead, 
ELLIPKE  uses  m  which  is  known  as  the  parameter.  The  modulus,  m,  is  related  to  the  parameter, 
k,  as  shown  in  Equation  3-14. 


—  m 

3-14 

The  induced  velocity  can  be  calculated  at  any  location  except  a  point  on  the  vortex  ring  (i.e.  the 
velocity  a  vortex  ring  induces  on  itself).  Future  versions  of  OpenProp  may  require  the  duct 
algorithm  to  calculate  the  velocity  a  vortex  ring  induces  on  itself.  Two  options  were  explored  to 
overcome  this  singularity  problem. 

First,  the  average  of  the  induced  velocity  of  selected  locations  on  a  sphere  surrounding  a  point  on 
the  vortex  ring  was  calculated.  This  was  done  for  consecutively  smaller  spheres.  The  result  did 
not  converge,  but  rather  grew  without  bound  as  the  sphere  size  was  reduced. 

The  second  attempt  to  overcome  the  singularity  problem  consisted  of  discretizing  the  vortex  ring 
into  vortex  filaments.  The  influence  of  each  vortex  filament  on  the  desired  point  was  summed. 
The  vortex  filament  that  contained  the  desired  point  was  not  included  in  this  summation  as  it  was 
assumed  that  the  point  was  on  this  filament,  and  it  is  well  know  that  there  is  no  influence  on  any 
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point  collinear  with  a  vortex  filament.  This  method  also  failed  to  eonverge.  As  the  discretation 
inereased,  so  did  the  resulting  indueed  veloeity. 

Table  1  gives  the  result  for  the  self-indueed  axial  veloeity  of  a  vortex  ring  of  unit  radius  and 
strength. 


Discretation 

Axial  Induced  Velocity 

10' 

0.2256 

W 

0.4072 

10' 

0.5904 

10" 

0.7737 

lo^ 

0.9569 

10^ 

1.1401 

10' 

1.3234 

I(? 

Not  Enough  Memory 

Table  1:  Vortex  Ring  Self-Induced  Axial  Velocity 


The  failure  to  overeome  this  singularity  did  not  affeet  the  duet  algorithm  used  in  OpenProp  v2.0 
beeause  the  self-infiuenee  of  a  ring  vortex  is  not  required.  If  it  is  required  in  a  future  version, 
the  error  ean  be  mitigated  by  increasing  the  discretation  of  the  duct  (i.e.  the  number  of  vortex 
rings  used  to  model  the  duct  is  increased). 


^  Real  vortices  do  not  have  this  singularity  problem  because  of  their  viscous  core  dissipation. 
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3.3  Circumferential  Mean  Velocity 


When  a  propulsor  includes  more  than  one  component,  it  becomes  necessary  to  calculate  the 
velocities  that  one  component  induces  on  another.  For  rotating  components  such  as  the 
propeller,  the  time-averaged  induced  velocities  are  used  and  are  equal  to  the  circumferential 
mean  velocities  calculated  in  the  rotating  reference  frame  of  the  component.  Formulas  for 
calculating  the  tangential,  axial,  and  radial  induced  velocities  induced  from  a  horseshoe  vortex 
are  presented  below.  The  formulas  are  derived  from  Coney  (1)  and  Hough  and  Ordway  (7),  and 
the  notation  most  closely  matches  Coney  (1). 


From  Kelvin’s  theorem.  Equation  3-15  gives  the  tangential  circumferential  mean  velocity,  ul, 
induced  on  a  control  point  at  radius  rc(m)  of  another  component  from  a  horseshoe  vortex  of 
strength  F with  lattice  points  at  radii  Vyiy)  —  1)  and  r^,(p).  Equation  3-16  defines  the  parameter  S 
which  directly  relates  to  whether  or  not  the  control  point  is  in  the  slipstream.  Xf  is  the  axial 
distance  from  the  horseshoe  vortex  lattice  point  to  the  control  point  with  positive  Xf  being  in  the 
downstream  (i.e.  positive  x-axis)  direction. 


ul  (m,  n) 


0, 

0, 

-zr 


^Znrfrn)  ’ 


S  >  0,  —  oo  <  Xf  <  oo 

S  <  0,  Xf  <  0 

S  <  0,  Xf  >  0 


3-15 


5  =  (r^(p  -  1)  -  rc(m))(r^(p)  -  rfm)) 

3-16 

That  is,  the  tangential  velocity  vanishes  everywhere  outside  of  the  slipstream  of  the  horseshoe 
vortex  and  is  proportional  to  r(p)/rc(m)  inside  the  slipstream.  The  tangential  circumferential 
mean  velocity  induced  by  both  the  bound  and  trailing  vorticity  can  be  found  from  the  above 
equations. 
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The  bound  vortices  on  a  radial  lifting  line  only  induce  tangential  circumferential  mean  velocities. 

The  axial  and  radial  circumferential  mean  velocities  induced  from  the  trailing  vortices  must  now 

be  calculated.  Hough  and  Ordway  (7)  used  Fourier  analysis  to  derive  formulas  for  the  induced 

velocities  in  terms  of  the  Heuman  Lambda  function  and  Legendre  functions  of  the  second  kind 

and  half  integer  order.  As  Coney  noted  in  (1),  these  can  be  thought  of  as  the  velocities  induced 

by  a  propeller  with  an  infinite  number  of  blades,  and  since  the  circumferential  mean  velocities 

are  the  average  of  the  sum  of  local  induced  velocities  along  a  circle,  Equations  3-17  and  3-18  can 

be  applied  to  calculate  the  axial  and  radial  circumferential  mean  velocities.  The  constant  C;  is 

defined  in  Equation  3-19  where  Q.i(_q)  are  the  Eengendre  functions  of  the  second  kind  and  half 

-2 

integer  order  and  AqC^,  k)  is  the  Henman’s  Eambda  function  with  the  amplitude,  (p,  and 
modulus,  K,  as  the  arguments.  The  arguments  for  the  Eegendre  and  Heuman  Eambda  functions 
are  given  in  Equations  3-20,  3-21,  and  3-22. 


Uaim,p)  =  r 


7rr^(p)tan()g^(p)) 


Cl 


Uy(rn,p)  =  f- 


■  Qi 


n  rc(jn)r^(p)tan{p^(p)) 


Xf  n 

n  +  ,  ==Qli(d  +  r,(m)  <  r^p) 

2^r^(m)r„(p)  ^  / 

Xf  n 

r  ,  '  =9_i(7)  - r^(m)  >  r^p) 

2^r,{m)r^(p)  2  2: 


3-17 


3-18 


3-19 


q  =  1  + 


2 

xj  +  (r^Cm)  -  r^Cp)) 

2rc(m)r^(p) 


3-20 
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q)  =  sin" 


Xf 


xj  +  (rc(m)  -  ry(p)y 


3-21 


K  = 


4rc(m)rv(p) 

2 

+  (rc(m)  +  r^Cp)) 


3-22 


As  previously  stated,  the  bound  vortieity  of  eaeh  horseshoe  vortex  only  induees  a  tangential 
veloeity.  Therefore,  equations  3-17  and  3-18  ean  be  applied  as  shown  in  Equations  3-23  and 
3-24  to  ealeulate  the  axial  and  radial  veloeities  indueed  from  a  horseshoe  vortex  that  is  used  in 
representing  a  radial  lifting  line. 


u*(m,p)  =  Uaim,p)  -  Ua(jn,p  +  1) 


3-23 


Ur(jn,p)  =  Ur(m,p)  —Uy(jn,p  +  1) 


3-24 


For  OpenProp,  the  above  algorithms  were  implemented  in  the  CMV.m  function  (Appendix  A). 
This  function  was  validated  with  Tables  1  and  2  in  (7).  The  test  case  used  the  representative 
blade  circulation  distribution  shown  in  Figure  3-4  and  assumed  that  the  helical  path  of  the 
trailing  vortex  system  was  determined  solely  by  the  incoming  free  stream  and  propeller  rotation 
(i.e.  the  advance  angle,  P,  was  used  instead  of  the  hydrodynamic  advance  angle.  Pi). 
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Figure  3-4:  CMV.m  Validation  Representative  Circulation 


CMV.m  required  two  Legendre  functions  and  the  Heuman  Lambda  function.  As  these  functions 
were  not  available  as  class  functions  in  MATLAB®,  they  were  written  and  validated  for  this 
thesis.  They  are  described  below. 

Q2half.m  (Appendix  B)  computes  the  Legendre  function  of  the  second  kind  and  positive  half 
order  of  the  argument  q  in  accordance  with  (8)  as  shown  in  Equation  3-25. 


3-25 

K(k)  denotes  the  complete  elliptic  integral  of  the  first  kind,  and  E(k)  denotes  the  complete  elliptic 
integral  of  the  second  kind  with  the  modulus,  k,  as  the  argument.  This  was  implemented  in 
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MATLAB®  using  ELLIPKE  with  the  parameter,  m,  as  the  argument  whieh  is  shown  in  Equation 
3-26. 


m  =  = - 

q  +  1 


3-26 


This  function  was  validated  with  Table  XIII  in  (9).  For  example: 

QlCl.S)  =  Q2half  (l.S)  =  0.39175 

2 

Qi(2.7)  =  Q2half  (2.7)  =  0.134035 

2 

Qi(8.4)  =  Q2half  (8.4)  =  0.0229646 

2 


Q2Mhalf.m  (Appendix  B)  computes  the  Legendre  function  of  the  second  kind  and  minus  half 
order  of  the  argument  q  in  accordance  with  (8)  as  shown  in  Equation  3-21 . 


3-27 

As  before,  K(k)  denotes  the  complete  elliptic  integral  of  the  first  kind  with  the  modulus,  k,  as  the 
argument.  This  was  implemented  in  MATLAB®  using  ELLIPKE  with  the  parameter,  m,  as  the 
argument  as  described  above  for  Qlhalf.m. 

This  function  was  validated  with  Table  XIII  in  (9).  For  example: 

(2_i(1.5)  =  Q2half  (l.S)  =  2.01891 

2 

<2_l(2.7)  =  Q2half  (2.7)  =  1.38958 

2 

<2_l(8.4)  =  Q2half  (8.4)  =  0.768523 

2 
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Heuman.m  (Appendix  B)  computes  Heuman’s  Lambda  function  of  the  arguments  (p  {amplitude) 
and  a  {modular  angle)  in  accordance  with  (8)  as  shown  in  Equation  3-28. 


2 

=  -{K-ic^)E{(p^  -a)  -  [K{a)  -  E(a)]F((p\|  -  a)} 


3-28 


a  =  sin  ^(k) 


3-29 


K(a)  denotes  the  complete  elliptic  integral  of  the  first  kind  with  the  modular  angle,  a,  as  the 
argument.  F((p\^  —  a)  denotes  the  incomplete  elliptic  integral  of  the  first  kind  with  the 
amplitude,  (p,  and  complementary  modular  angle,  |  —  a,  as  the  arguments.  E((p\^  —  a)  denotes 
the  incomplete  elliptic  integral  of  the  second  kind,  with  the  amplitude,  (p,  and  complementary 
modular  angle,  |  —  a,  as  the  arguments. 

K(a)  was  implemented  in  MATLAB®  using  ELLIPKE  with  the  parameter,  m,  as  the  argument 
as  shown  in  Equation  3-30. 


m  =  —  sin  (a)^ 

3-30 

F((p\|  —  a)  was  implemented  in  MATLAB®  using  the  imbedded  Maple  function  F/h)!7hcF  with 
the  sine  of  the  amplitude,  sin((p),  and  the  parameter,  k,  as  the  arguments.  Since  the 
complementary  modular  angle  is  used  for  the  incomplete  elliptic  integral  in  this  case,  the 
parameter  is  defined  as  shown  in  Equation  3-31. 

k  —  sin  —  aj 

3-31 
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E((p\|  —  a)  was  implemented  similarly  in  MATLAB®  using  the  imbedded  Maple  function 
EllipticE. 


Eleuman.m  was  validated  with  (8). 

Figure  3-5  was  generated  using  Eleuman.m  and  agrees  with  Figure  17.10  in  (8). 


Figure  3-5:  Heuman’s  Lambda  Function 
The  sample  calculations'^  given  below  agree  with  Table  17.8  in  (8). 


ylo(5°\10°)  =  HeumaniS°\10°)  =  0.086495 
ylo(45°\60°)  =  Heuman(4S°\60°)  =  0.569122 
ylo(75°\40°)  =  Heuman(7S°\40°)  =  0.906056 


*  When  executing  Heuman.m  in  MATLAB,  cp  and  a  are  entered  in  radians  vice  degrees. 
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3.4  Circulation  Optimization 

The  goal  of  OpenProp’s  optimization  routine  is  to  calculate  the  radial  distribution  of  circulation 
on  the  lifting  line  that  minimizes  torque  for  a  given  thrust.  Also  specified  are  the  propeller 
diameter,  number  of  blades,  advance  coefficient  (Js),  and  inflow  velocity  profile.  OpenProp  vl.O 
uses  the  Lerbs  criterion  where  tan()ffj)  is  obtained  from  tan()ff)  in  terms  of  an  unknown 
multiplicative  factor  (10).  The  optimization  routine  initially  estimates  the  hydrodynamic  pitch 
angle  (Pi)  based  on  the  undisturbed  flow  angle  (P)  and  the  efficiency  of  the  actuator  disk.  The 
system  of  equations  represented  by  Equation  3-32  is  then  solved  to  obtain  the  optimum 
circulation: 


M 

^  [uain,m)  -  Utin,m)  tan(p(n))]  =  Va(n) 

m=l 


/tan(pi(n))  \ 
Vtan(p(n))  ) 


n  =  1,  ...M 


3-32 


Using  the  circulation  distribution,  the  vortex  induced  velocities  on  each  panel  of  the  lifting  line 
are  then  solved.  This  allows  the  forces  to  be  calculated,  and  the  resulting  thrust  is  compared  with 
the  desired  thrust.  The  hydrodynamic  pitch  angle  is  then  iteratively  adjusted,  and  the  process  is 
repeated  until  the  desired  thrust  is  achieved  (2). 


This  thesis  integrated  the  duct  model  into  the  Lerbs-based  optimization  routine  as  follows: 

The  influence  functions  included  the  effects  of  the  duct  image  horseshoes. 

The  inflow  velocities  were  modified  to  include  the  effect  of  the  duct  vortex  rings. 

The  duct  circulation  was  an  entered  value.  The  circumferential  mean  velocity 
induced  by  the  propeller  on  the  duct  was  added  to  the  inflow,  and  the  thrust  produced 
by  the  duct  was  calculated  using  the  Kutta-Joukowski  law. 

The  duct  thrust  was  subtracted  from  the  desired  thrust. 


Using  PEL,  several  test  cases  were  run  for  validation.  In  each  case,  the  results  indicated  that  the 
existing  Eerbs-based  optimization  routine  in  OpenProp  did  not  converge  to  the  same  result  as 
PEE  which  uses  a  calculus  of  variations  optimization  routine.  The  results  for  the  t  =  0.8  case 
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for  a  5-bladed,  10  foot  diameter  propeller  with  Ct=1.2  and  Js=0.60  is  presented  below.  Viscous 
forces  were  ignored. 


Radius  (r/R) 


Figure  3-6:  OpenProp  t=0.80  test  case  using  the  Lerbs-based  optimization  routine 


The  results  showed  that  the  Lerbs-based  optimization  routine  in  its  current  form  was  not 
extendable  to  a  ducted  propeller.  Upon  further  consideration,  this  was  logical  given  that  the 
routine  relied  on  the  initial  undisturbed  flow  angle,  P,  at  the  propeller  plane  but  did  include  the 
effect  of  the  duct  in  this  initial  determination  of  p.  It  may  be  possible  to  integrate  a  P  iteration 
loop  into  the  Lerbs-based  routine,  but  that  option  was  not  pursued  for  this  thesis. 


Following  the  unsuccessful  attempt  with  the  Lerbs-based  optimization  routine,  the  variational 
optimization  routine  presented  in  (1)  was  integrated  into  OpenProp.  This  option  was  chosen  for 
two  reasons:  1)  it  is  the  proven  routine  used  in  PLL  and  2)  it  is  a  more  general  optimizer  that  can 
be  used  for  many  different  propulsor  configurations.  In  this  routine  the  wake  geometry  is  frozen 
while  the  optimum  circulation  distribution  is  calculated  such  that  torque. 
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M 

Q  =  pz  1  [14  (m)  +  u*  (m)]r(m)r(m)Ar, 

m=l 

3-33 

is  minimized,  subject  to  the  constraint  that  the  thrust, 

M 

T  =  pZ  1  [Vf-(m)  +  a)r(m)  -I-  ul(jn)]T(m)Ar, 

m=l 

3-34 

has  a  prescribed  value,  T^.. 


The  auxiliary  function  H  =  Q  +  A(T  —  is  formed,  and  its  partial  derivatives  with  respect  to 
the  unknown  variables  (  r(i)'s  and  X  )  are  set  to  zero  as  shown  in  equations  3-35  and  3-36 
below.  The  Lagrange  multiplier,  1,  is  an  additional  unknown  variable  and  must  be  solved  for 
along  with  the  discrete  circulation  strengths,  the  r(i)'s. 


dH 


for  i  =  1..M 


3-35 

3-36 

Equations  3-35  and  3-36  form  a  nonlinear  system  of  M  +  1  equations  with  M  unknown  values  of 
circulation  and  an  unknown  Lagrange  multiplier.  By  assuming  that  the  Lagrange  multiplier  is 
known  where  it  forms  quadratic  terms  with  the  circulation  and  that  the  tangential  induced 
velocity,  ul(ra),  is  known,  the  solution  to  the  nonlinear  system  of  equations  can  be  found  by 
iteratively  solving  the  following  linear  system  of  equations: 
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dH 

Wio 


0  =  Vaii)r(i)Ar 


M 


+ 


^  [r(m)ua(i,  m)r(m)Ar  +  r(m)ua(m,  i)r(i)Ar] 


m=l 


+A[Vf-(i)  +  a)r(i)]Ar 


M 


+X  ^  [r(m)ut(i,  m)Ar  +  r(m)ut(m,  i)Ar] 


m=l 


for  i  —  1  ...M 


3-37 


M 

Tr^pZ  1  [Vf-(rn)  +  a)r(m)  +  ul(m)]T(m)Ar, 

m=l 

3-38 

For  each  iteration,  the  frozen  Lagrange  multiplier,  X,  in  equation  3-37  and  the  tangential  induced 
velocity  in  equation  3-38  take  on  values  from  the  previous  values.  Furthermore,  Coney 
determined  that  initially  setting  the  induced  velocities  equal  to  zero  and  the  Lagrange  multiplier 
equal  to  -1  are  suitable  initial  estimates  of  these  quantities. 

After  the  optimum  circulation  distribution  is  found,  the  wake  is  aligned  with  the  velocities 
induced  by  that  circulation  distribution.  The  aligned  wake  is  now  frozen  and  the  above  process 
is  repeated  in  order  to  determine  a  new  optimum  circulation  distribution.  By  using  this  double 
iterative  approach,  velocities  and  forces  consistent  with  moderately  loaded  lifting  line  theory  can 
be  obtained. 

The  above  procedure  was  adapted  to  handle  a  ducted  propeller.  Additionally,  an  estimate  of  the 
effects  of  viscous  drag  was  added.  The  following  modifications  were  made: 

The  effects  of  the  duct  image  horseshoes  were  added  to  the  influence  functions  in 
equation  3-37  as  shown  in  equations  3-6  and  3-7. 

The  effects  of  the  duct  ring  vortices  were  included  in  calculating  the  induced 
velocities  on  the  lifting  line. 
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The  NACA  a=0.8  meanline  eireulation  distribution  for  the  duct  ring  vortex  system 
was  scaled  so  that  the  duct  produces  the  thrust  specified  by  the  thrust  ratio. 

The  duct  thrust,  T^,  calculated  as  shown  in  equation  3-9,  was  subtracted  from  the 
required  thrust  in  equation  3-38.  Additionally,  an  estimate  of  the  propeller’s  viscous 
drag,  Tyiscous  >  was  added  to  the  required  thrust.  The  modified  required  thrust  is 
given  in  Equation  3-39. 

M 

Tr=pZ^  [VtOn)  +  mr(m)  +  u^(m)]r(m)Ar  +  Tyiscous  “  Ta 

m=l 

3-39 


The  propeller’s  viscous  drag  was  estimated  using  the  two-dimensional  airfoil  sectional  drag 
coefficients,  CD^(jn),  and  the  section  chord  lengths,  c(m)  as  shown  in  Equation  3-40. 

1  ^ 

T'viscous  ~  1 

m=l 

3-40 


In  OpenProp  v2,  the  variational  optimization  routine  described  above  was  integrated  as  a 
MATEAB®  function  named  Coney.m^  (Appendix  C).  Within  Coney.m,  the  ductVort.m  and 
ductThrust.m  functions  discussed  above  were  used  to  perform  specific  duct-related  calculations. 


^  Brenden  Epps  wrote  the  initial  version  of  Coney.m.  The  version  of  Coney.m  implemented  in  OpenProp  v2  was  a 
modified  and  expanded  version  of  his  eode. 
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4.  Implementation  and  Validation 

In  order  to  enable  OpenProp  to  design  a  dueted  propeller,  the  variational  optimization  routine 
{Coney. m)  and  associated  MATLAB®  functions  {ductVort.m,  ductThrust.m,  vRing.m,  CMV.m, 
Qlhalf.m,  QlMhalf.m,  and  Heuman.m)  discussed  above  were  integrated  into  the  OpenProp  code. 
The  graphical  user  interfaces  (GUIs)  were  updated  to  include  the  required  parameters  for  a 
ducted  propeller  as  shown  in  Figure  4-1. 


Figure  4-1:  OpenProp  v2  single  propeller  design  GUI  with  duct  parameters 


The  required  duct  parameters  are  located  in  the  upper-right  corner  and  are  shown  enlarged  in 
Figure  4-2. 


I~  Hub  Image  Flag  (Check  for  YES) 
R  Ducted  propeller  (Check  for  YES) 


nr 

nr 

noo8 


Thrust  Ratio 

Duct  Diameter/Prop  Diameter 
Duct  Section  Drag  Coefficient 


Figure  4-2:  OpenProp  v2  required  duct  parameters 


If  the  user  desires  to  design  a  propeller  operating  within  a  duct,  the  Ducted  Propeller  check  box 
is  selected  and  three  parameters  are  required: 
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Thrust  Ratio:  as  given  in  equation  3-2,  the  thrust  ratio  specifies  the  portion  of  the 
desired  thrust  that  the  propeller  will  produce.  Typical  values  of  thrust  ratio  are  0.7  to 
1 .3.  Values  less  than  one  result  in  an  accelerating  duct  that  produces  a  positive  thrust 
while  values  greater  than  one  result  in  a  decelerating  duct  that  produces  a  negative 
thrust  which  must  be  overcome  by  the  propeller. 

Duct  Diameter/Prop  Diameter:  determines  the  size  of  the  duct  in  terms  of  the 
propeller  diameter.  This  value  must  be  equal  to  or  greater  than  one  as  a  duct  must  at 
least  be  as  large  as  the  propeller  which  it  surrounds.  This  parameter  determines  the 
gap  between  the  duct  and  the  propeller  tip  which  has  a  major  influence  on  the 
optimum  circulation  distribution.  For  this  thesis,  the  Duct  Diameter/Prop  Diameter 
parameter  was  set  to  one  and  disabled  since  only  the  zero  gap  case  was  considered. 
Duct  Section  Drag  Coefficient:  specifies  the  two-dimensional  airfoil  sectional  drag 
coefficient  which  is  used  to  estimate  the  duct’s  viscous  drag.  This  parameter  must  be 
equal  to  or  greater  than  zero.  A  value  of  zero  implies  the  inviscid  case  where  viscous 
effects  are  neglected.  A  typical  value  for  the  viscous  case  is  0.008. 

OpenProp  v2  assumes  the  following: 

The  duct  surrounds  a  radial  lifting  line. 

The  duct  chord  length,  q,  is  equal  to  the  propeller  radius,  R. 

The  duct  is  positioned  such  that  the  mid-chord  is  located  at  the  lifting  line. 

The  circulation  distribution  of  the  duct  vortex  ring  system  represents  a  NACA  a=0.8 
meanline. 

In  order  to  demonstrate  the  capability  of  OpenProp  v2  and  provide  validation  via  a  PLL 
comparison,  several  test  cases  were  completed  using  a  5-bladed  ducted  propeller  with  optimum 
circulation  distribution  operating  at  Js  =  0.60  and  Cj  =  1.20.  The  duct  had  zero  thickness,  and 
there  was  zero  gap  between  the  duct  and  the  propeller.  For  viscid  runs,  a  sectional  drag 
coefficient  of  0.008  was  used  for  both  the  duct  and  the  propeller.  Appendix  D  contains  all  of  the 
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run  settings  for  both  OpenProp  v2  and  PLL^.  Figure  4-3  shows  a  sample  rendering  of  the  ducted 
propellers  designed  in  the  test  cases. 


Figure  4-3:  Sample  rendering  of  ducted  propeller  produced  in  the  test  cases 


First,  to  ensure  the  various  functions  and  routines  in  the  optimization  routine  were  working 
correctly,  the  circulation  from  a  t  =  0.8  PLL  run  was  fed  into  OpenProp  v2.  Figure  4-4  and 
Figure  4-5  show  given  the  same  circulation  distribution,  OpenProp  v2  calculated  the  same 
induced  velocities  as  PLL.  This  included  the  total  induced  velocities  on  the  propeller  control 
points  (UA*/Vs  and  UT*/Vs),  the  axial  velocities  induced  by  the  duct  rings  on  the  propeller 
control  points,  and  the  velocities  on  the  duct  rings  by  the  propeller  lifting  line. 


®  All  runs  for  both  PLL  and  OpenProp  v2  used  10  vortex  panels.  Neither  PLL  or  OpenProp  v2  would  run 
sueeessfully  with  greater  than  20  vortex  panels.  PLL  would  erash  for  an  unknown  reason,  and  OpenProp  v2 
experieneed  an  error  in  the  Heuman.m  funetion. 
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UAWS 


Tau=0.8  (inviscid)  Tau=0.8  (inviscid) 


Tau=0.8  (inviscid) 


Tau=0.8  (inviscid) 


Figure  4-4:  OpenProp  v2  algorithm  propeller  results  usiug  PLL  circulatiou 
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Tau=0.8  (inviscid) 


Tau=0.8  (inviscid) 


Tau=0.8  (inviscid)  Tau=0.8  (inviscid) 


Figure  4-5:  OpenProp  v2  algorithm  duct  results  usiug  PLL  circulatiou 


Figure  4-6  gives  an  efficiency  versus  thrust  ratio  comparison  for  PLL  and  OpenProp  v2.  The 
ideal  efficiency  as  calculated  by  equation  3-1  using  actuator  disk  theory  is  also  shown.  For  this 
figure,  PLL  and  OpenProp  v2  each  ran  independently  (i.e.  OpenProp  v2  calculated  its  own 
optimum  circulation  distribution  instead  of  using  PLL’s  as  was  the  case  in  Figure  4-4  and  Figure 
4-5). 
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Figure  4-6:  Efficiency  versus  thrust  ratio  (t)  comparison  between  OpenProp  v2  and  PEL 


The  OpenProp  v2  results  did  not  match  PLL  exactly,  but  they  are  close  (within  1%  for  t  <  1.0) 
and  follow  the  same  trend  as  PLL.  As  with  PLL,  the  maximum  efficiency  for  the  ducted 
propeller  occurred  at  a  thrust  ratio  of  approximately  t  =  0.9.  The  difference  between  OpenProp 
v2  and  PLL  increased  above  t  =  1.0  and  was  approximately  2%  at  t  =  1.2.^  The  reason  for  the 
difference  is  explained  by  the  optimum  circulation  distribution  calculated  by  OpenProp  v2. 

The  following  group  of  figures  show  comparisons  between  OpenProp  v2  and  PLL  for  thrust 
ratios  of  0.8,  1.0,  and  1.2.  Viscosity  is  neglected.  For  t  =  0.8,  three  figures  are  shown:  Figure 
4-7  gives  the  standard  OpenProp  graphical  output  and  Figure  4-8  and  Figure  4-9  show  the 
comparison  between  OpenProp  v2  and  PLL.  For  t  =  1.0  and  t  =  1.2,  only  the  quad-chart  with 
the  circulation  comparison  is  shown. 


^  OpenProp  v2  did  not  converge  for  thrust  ratios  greater  than  1.2. 
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Non-Dimensional  Circulation 


J=0.600;  Ct=1 .200;  Cq=0.295;  Kt=0.170:  Kq=0.021 ;  ri=0.776:  x=0.800 


1.2 
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Figure  4-7:  OpenProp  v2  graphical  output  for  test  case  with  t  =  0. 8 


Tau=0.8  (inviscid)  Tau=0.8  (inviscid) 


Tau=0.8  (inviscid)  Tau=0.8  (inviscid) 


Figure  4-8:  OpenProp  v2  comparison  with  PLL  for  t  =  0. 8 
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Tau=0.8  (inviscid) 


Q 

O 


Tau=0.8  (inviscid) 


Tau=0.8  (inviscid) 


Figure  4-9:  OpenProp  v2  comparison  with  PLL  for  t  =  0. 8  (duct  ring  velocities) 

Figure  4-8  shows  that  the  optimum  eireulation  distribution  obtained  from  OpenProp  v2  is 
slightly  different  than  PLL’s  solution  for  the  accelerating  duct  case.  The  main  difference  was 
that  with  OpenProp  v2  the  circulation  reached  a  maximum  value  at  approximately  r/R  =  0.7  and 
decreased  slightly  at  the  tip. 


Figure  4-10  shows  the  comparison  between  OpenProp  v2  and  PLL  for  t  =  1.0  (neutral  duct). 
The  results  match  very  well. 
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Tau=1 .0  (inviscid) 


Tau=1 .0  (inviscid) 


Tau=1 .0  (inviscid)  Tau=1 .0  (inviscid) 


Figure  4-10:  OpenProp  v2  comparison  with  PLL  for  t  =  1. 0 


Figure  4-11  shows  the  comparison  between  OpenProp  v2  and  PLL  for  the  decelerating  duct  case 
of  T  =  1.2.  As  with  the  accelerating  duct  case,  the  optimum  circulation  distribution  obtained 
from  OpenProp  v2  is  slightly  different  than  PLL’s  solution.  For  the  decelerating  duct,  OpenProp 
v2’s  optimum  circulation  distribution  has  an  inflection  point  at  approximately  r/R  =  0.7  and  the 
tip  circulation  is  slightly  higher  than  the  PLL  solution. 
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UAWS 


Tau=1 .2  (Inviscid) 


Tau=1 .2  (inviscid) 


Tau=1 .2  (inviscid)  Tau=1 .2  (inviscid) 


Figure  4-11:  OpenProp  v2  comparison  with  PLL  for  t  =  1. 2 
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5.  Conclusions  and  Recommendations 


5.7  Conclusions 

This  thesis  successfully  extended  OpenProp’s  capability  such  that  it  can  now  design  a  propeller 
operating  inside  a  duct  with  no  gap  for  thrust  ratios  between  0.7  and  1.2.  The  results  confirmed 
that  maximum  efficiency  is  obtained  with  a  thrust  ratio  of  approximately  0.9. 

The  variational  optimization  routine  used  in  OpenProp  v2  was  validated  with  output  from  PLL. 
All  induced  velocity  calculations  matched  PLL,  and  the  optimum  circulation  distribution 
matched  PLL  for  the  neutral  duct  case.  However,  the  optimum  circulation  distribution  obtained 
from  OpenProp  v2  for  both  the  accelerating  and  decelerating  duct  cases  varied  slightly  but 
distinctively  from  PLL.  With  PLL,  the  circulation  distribution  always  reached  a  maximum  at  the 
tip.  This  was  not  the  case  with  OpenProp  v2.  Only  for  the  neutral  duct  case  was  this  true  for 
OpenProp  v2.  For  the  accelerating  duct,  the  circulation  peaked  at  approximately  r/R  =  0.7  and 
the  tip  circulation  was  lower  than  PLL’s  tip  circulation.  For  the  decelerating  duct,  the  circulation 
had  an  inflection  point  at  approximately  r/R  =  0.7  and  the  tip  circulation  was  higher  than  PLL’s. 

The  specific  reason  for  the  differing  optimum  circulation  distributions  was  not  discovered. 
However,  the  author  did  verify  that  the  variational  optimization  routine  implemented  in 
OpenProp  v2  was  a  faithful  representation  of  the  routine  presented  by  Coney  in  (1).  It  was 
assumed  that  PLL  used  this  variational  optimization  routine  as  well,  but  it  is  possible  that  PLL 
added  additional  constraints  that  were  not  discussed  in  (1). 
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5.2  Recommendations  for  further  work 


This  thesis  only  examined  a  duct  with  no  thickness.  Source  rings  could  be  integrated  into  the 
algorithm  to  model  duct  thickness. 

A  tip  gap  model  could  be  added  to  represent  the  flow  between  the  propeller  and  duct  when  a  gap 
exists.  This  would  allow  the  design  of  ducted  propellers  with  gaps  greater  than  zero. 

This  thesis  did  not  analyze  the  flow  around  the  duct  nor  did  it  attempt  to  define  the  true  duct 
orientation.  By  analyzing  the  flow  around  the  duct  three  important  objectives  could  be  obtained. 
First,  the  designer  could  ensure  that  flow  separation  does  not  occur  on  the  duct.  This  is  critical 
because  if  separation  occurs,  drag  will  increase  dramatically.  Second,  understanding  the  flow 
characteristics  in  the  gap  is  essential  to  analyzing  the  performance  of  the  entire  system  under 
various  loading  conditions.  Third,  the  flow  streamlines  would  outline  the  shape  of  the  duct  and 
reveal  the  duct  angle  of  attack  (i.e.  orientation).  Coupling  OpenProp  with  a  computational  fluid 
dynamics  code  would  be  the  ultimate  goal  to  properly  analyze  the  flow. 
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Appendix  A.  Duct  Theory  MATLAB®  Code 


A.l  ductVortm 

%Discrete  representation  of  vorticity  (circulation)  on  vortex  rings 
%that  represent  duct  as  a  NACA  a=0 . 8  meanline. 


%Calculates  duct  vortex  ring  influence  (UADUCT) 
%on  prop  lifting  line  Ctrl  pts 


%Variables : 

%  %  R  [m]  : 

%  %  rDuct 

g,  q, 
o  o 

%  %  Mp: 

O,  o,  pp  , 

o  o  rxL/  • 


[m]  : 


propeller  radius 

duct  radius  (formerly  vrRad) ,  vortex  ring  radius  at 
xDuct 

#  of  control  points  (N=M,  M:  #  of  vortex  rings) 
radius  of  Ctrl  pts  on  lifting  line 


%Note:  duct  chord  =  R 


%Returns : 

q, 

o 

g, 

o 

g, 

o 

q, 

o 

g, 

o 


vRingLoc  [m] : 

dVort  [m''2/s]  : 
UADUCT  [m/s] : 


location  of  each  vortex  ring  (x,y,z  vector), 
(formerly  (vortex) 

circulation  distribution  of  each  vortex  ring 
duct  vortex  ring  axial  influence  on  prop 
lifting  line  control  pts 


%close  all;clear  all;clc; 


function  [vRingLoc, dVort, UADUCT]  =  ductVort (R, rDuct, Mp, RC) 

%setup  vRing  spacing 
M=Mp+2; 

if  rem(M,2)~=0  %ensures  Mp  is  even 

M=M+1; 

end 

dS=l/M;  %spacing  between  vRings  (non-dim  with  R) 

hdS=0.5*dS;  %half  of  dS 


%computes  the  circulation  on  vortex  rings  which  each  represent  the 
%vorticity  on  a  piece  of  NACA  a=0.8  mean  line  located  at  position  XvRing 
%(L.E.=0.0,  T.E.=1.0)  and  is  of  length  dS . 

XvRing=zeros (1,M) ; dVort=XvRing; 

for  n=l : M 

XvRing (n) = (n-1) *dS+hdS; 

X2  =  XvRing (n)  +  hdS; 

XI  =  XvRing (n)  -  hdS; 
if  X2  <=  0.8 

dVort (n)  =  dS/0.9; 
elseif  XI  >=  0.8 

Y1  =  1.0  -  (XI  -  0.8) /0.2; 
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Y2  =  1.0  -  (X2  -  0.8) /0. 2; 
dVort(n)  =  dS*0.5*(Yl  +  Y2)/0.9; 

else 

Y2  =  1.0  -  (X2  -  0.8) /0. 2; 

FRONT  =  0.8  -  XI; 

BACK  =  0.5* (1.0  +  Y2)*(X2  -  0.8); 
dVort(n)  =  (FRONT  +  BACK) /O. 9; 

end 

end 

LED=- (M/2 ) *dS ;  %location  of  leading  vRing  on  duct 

XvRing=XvRing+LED; 

vRingLoc=zeros (3,M)  ; 

vRingLoc ( 1 , : ) =XvRing; 

vRingLoc(2,  ;)=rDuct/R; 

%  %  plot (XvRing, dVort, ' * ' , XvRing, dVortl , ' + ' ) 

%Calc  duct  vortex  ring  influence  (UADUCT)  on  prop  lifting  line  Ctrl  pts 
%Note:  No  tangential  influence 

%Note:  Radial  influence  does  not  create  a  force  on  radial  lifting  line 
UADUCT=zeros (Mp, 1 ) ;  %axial  influence 

URDUCT=UADUCT;UTDUCT=UADUCT;  %radial  and  tantential  influence 

for  n=l:Mp  %cycle  thru  all  Ctrl  pts  on  lifting  line 

P= [ 0 ; RC (n) ; 0 ] ;  %3D  coord  for  Ctrl  pt 

for  m=l:M  %cycle  thru  all  vortex  rings  on  duct 

UD  =  vRing (vRingLoc ( 1 , m) , vRingLoc (2 , m) , P, dVort (m) ) ; 

UADUCT (n) =UADUCT (n) +UD ( 1 )  ; 

URDUCT (n) =URDUCT (n) +UD (2) ; 

UTDUCT (n) =UTDUCT (n) +UD (3) ; 

end 

end 

UADUCT=UADUCT*2 *pi ;  %2*pi  needed  for  non-dimensional  circulation  (G) 

URDUCT=URDUCT*2*pi; 

UTDUCT=URDUCT*2*pi; 
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A.2  ductThrustm 


%Calculates  total  duct  thrust  coefficient  and  associated  parameters 

%This  version  handles  one  propeller.  Modifications  required  if  additional 
%rotors  or  stators  desired. 


g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 


Variables : 

%  vRingLoc: 

%  dVort: 

%  dCirc: 

%  rDuct: 

%  CDd: 

%  U  [m/s] : 

%  rho  [kg/m3] : 
%  RV: 

%  G: 

%  TanBI: 

%  Z: 

%  R: 


location  of  duct  vortex  rings 
circulation  distribution  on  duct  rings 
strength  of  duct  vorticity 
radius  of  duct 

coefficient  of  drag  for  the  duct 
x-dir  inflow  velocity  magnitude 
density  of  fluid 

radius  of  trailing  helical  vortices  on  lifting  line 
non-dim  circulation  for  panels  on  lifting  line 
tangent  of  betal  for  RV  points 
#  of  blades 
radius  of  propeller 


%Returns :  CTD: 

%  dCirc 
%  UAdVS 
%  URdVS 


total  duct  thrust  coeff  (viscous  drag  included) 
new  dCirc  scaled  to  provide  desired  duct  thrust 
induced  axl  velocity  on  duct  ring  from  lifting  line 
induced  rad  velocity  on  duct  ring  from  lifting  line 


function  [CTD, dCirc, UAdVS, URdVS]  =  ductThrust (vRingLoc, dVort, dCirc, .. . 

rDuct,  CDd,  U,  rho,  RV,  G,  TanBI,  Z,  R,  CTDDES) 


M=length (dVort) ; 
k=  [  0  0  1  ]  ; 

Uvec= [U; 0 ; 0 ] ; 


%#  of  duct  vortex  rings 
%unit  vector  in  Z  direction 
%free  stream  velocity  vector 


%Velocity  vector  at  each  vortex  (Vvortex) 
Vvortex=zeros (3,M) ; 
VvortexInduced=Vvortex; 

Lvortex=Vvortex; 


for  m=l;M  %cycles  thru  all  vortex  rings 

Vvortexinduced ( : , m) =VvortexInduced ( : , m)  +  CMVl (vRingLoc ( 1 , m) , ... 

vRingLoc (2, m) , RV, G, TanBI, Z) ;  %radial  lifting 
Vvortex (:,m)  =VvortexInduced ( : , m)  +  Uvec; 

Lvortex ( : , m)  =rho*cross (Vvortex ( : , m) , k*dCirc*dVort (m) ) ; 

%lift  on  a  vortex  ring 


end 


line 


URdVS=VvortexInduced (2 , : ) ;  %induced  rad  vel  on  duct  ring  from  prop 

%Note:  CMVl.m  values  for  axial  and  tangential  induced  velocities  on  the  duct 
%rings  don't  match  PLL  results.  Tangential  velocities  are  not  needed,  so  that 
%discrepancy  is  not  resolved.  Axial  velocities  only  match  for  -x  locations. 
%However,  from  PLL,  +x  are  a  negative  mirror  of  -x  values  so  I  have  adjusted 
%accordingly  so  that  axial  tangential  velocities  can  be  used  to  calculate  the 
%duct  viscous  drag. 
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%Adjust  axial  CMV  so  that  results  matches  PLL 
ml  =  M/2;  %last  -x  vortex 

for  m=l :M/ 2 

Vvortexinduced ( 1 , ml+m) =-VvortexInduced ( 1 , ml-m+1 ) ; 

end 

Vvortex ( 1 , : ) =VvortexInduced ( 1 , : )  +  Uvec(l); 

UAdVS=VvortexInduced ( 1 , : ) ;  %induced  axl  velocity  on  duct  ring  from  prop 


%Duct  thrust  (CONEY  P.  77,  EQN  3.28.) 

Lift=sum (Lvortex,  2 )  ; 

dThrust=-2*pi*rDuct*Lift (1) ;  %thrust  (-x  direction)  for  duct  [N]  :  '  ) 

%positive  implies  thrust  to  the  ship 


%Viscous  drag  for  duct  (Drag  =  0 . 5*rho*V''2*Chord*CDd  *  2*pi*rDuct) 
delS=abs (vRingLoc ( 1 , 1 ) -vRingLoc ( 1 , 2 ) ) ;  %vortex  spacing 

%linear  spacing  assumed 

dDrag  =  0; 
for  m=l :M 

dDrag  =  dDrag  +  Vvortex  ( 1 ,  m)  2  ; 


end 

dDrag=0 . 5*rho*dDrag*delS*CDd*2*pi*rDuct;  %R  =  duct  chord  length 

CTDdrag  =  4*dDrag  / (R*2*pi)/ (rho*R);  %normalized  duct  "drag"  CT 

CTDthrust  =  4*dThrust  / (rho*R);  %normalized  duct  "thrust"  CT 


dThrustTot=dThrust-dDrag/ (R*2*pi) ; 


%total  thrust  for  the  duct 
%(R*2*pi)  required  to  make 
%dDrag  dimensions  match 
%dThrust 


CTD  =  4*dThrustTot  /  (rho*R); 


%CT  for  duct 


%scale  duct  circulation  so  that  duct  provides  required  thrust 
if  dCirc~=0 

dCirc  =  dCirc/CTDthrust* (CTDDES+CTDdrag) ; 

end 
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A. 3  vRing.m 


%Returns  velocity  vector  induced  by  vortex  ring  of  input  strength 
%gamma  at  a  point  (p)  in  space  given  the  x-axis  location  (vrX)  and 
%radius  (vrRad)  of  the  vortex  ring: 

%Axis  of  the  vortex  ring  is  in  the  direction  of  the  x-axis. 


%Variables 
%gamma : 
%vrX: 
%vrRad : 

%p 


vortex  ring  strength 

x-axis  location  of  vortex  ring 

radius  of  vortex  ring 

(Px,Py,Pz)  point  at  which  velocity  is  induced 


%Returns :  Vp : 


velocity  at  point  P 


%Ref:  Kuchemann  and  Weber,  Aerodynamics  of  Propulsion  p  305. 


function  [Vp]  =  vRing (vrX, vrRad, p, gamma) 

if  vrRad  ==  0  %stops  function  if  vrRad  =  0 

Vp= [ 0 ;  0 ;  0 ] ; 
return 

end 

if  vrRad  <  0  %stops  function  if  vrRad  <  0 

Vp=[NaN;  NaN;  NaN] ; 
return 

end 


Px=p  (1)  ;  Py=p  (2)  ;  Pz=p  (3)  ; 

if  Pz==0 

if  Py<0 

thetaP=-pi/2 ; 

else 

thetaP=pi/2 ; 

end 

else 

if  Pz>0 

thetaP=atan (Py/Pz) ; 

else 

thetaP=atan (Py/Pz) +pi; 

end 

end 

Prad=sqrt (Pz^2  +  Py^2); 

if  vrX==p(l)  &  vrRad==Prad 
Vp=[NaN;  NaN;  NaN]  ; 
return 

end 

x= (Px-vrX) /vrRad; 
r=Prad/ vrRad; 


cylindrical  coord  angle  for  P 

logic  for  atan  ambiguity 

cylindrical  coord  radius  for  P 
stops  function  if  P  on  vortex  ring 

x/r'  from  Kuchemann 
r/r'  from  Kuchemann 
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%Elliptic  integral  method  (Kuchemann  p.  305) 

%uses  parameter  k  where  k^2  =  m  for  elliptic  integrals 

k=sqrt (4*r/ (x^2+ (r+1) ^2)); 

[K,  E]  =ellipke  (k''2)  ; 

Vx=gamma/  (2*pi*vrRad)  / sqrt  (x''2  +  (r+1)  "'2)  *  (K-  (1+2*  (r-1)  /  (x^2  +  (r-1)  "'2)  )  *E) 

if  r==0 
Vr=0; 

else 

Vr=gamma/ (2*pi*vrRad) * (-x) /r/ sqrt (x^2+ (r+1) ^2) * (K- (l+2*r/ . . . 
(x"2+(r-l) "2) ) *E) ; 

end 

Vy=Vr*sin (thetaP)  ; 

Vz=Vr*cos (thetaP) ; 

Vp=[Vx;  Vy;  Vz]; 
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A.4  vpfDuctm 


%Velocity  Prediction  Function  for  a  Duct 

%Returns  velocity  (Vp)  at  pt  P  due  to  a  set  of  vortex  rings, 
%source  rings  (for  thickness),  and  free  stream 


%Variables : 


P: 

vortex : 
gamma : 

S: 

Uvec : 


point  at  which  velocity  is  desired  (column  vector) 

matrix  of  vortex  and  source  locations 

matrix  of  vortex  ring  strengths 

matrix  of  source  ring  strengths 

free  stream  velocity  vector 


%Returns :  Vp : 


velocity  at  point  P 


function  [Vp]  =  vpfDuct (P, vortex, gamma, S, Uvec) 


Vp  =  [  0 ;  0 ;  0  ]  ; 

for  m=l : size (vortex, 2 )  %cycle  thru  all  vorticies 

Rvp=P-vortex ( : , m) ;  %vector  from  vortex  to  P 

if  norm (Rvp) <1 Oe-5  %skips  vortex  if  P  is  on  vortex 

else 

Jp=Rvp  /  (2  *pi*norm  (Rvp)  "^2  )  ; 

Vp=Vp  +  vRing (vortex ( 1 , m) , vortex (2 , m) , P, gamma (m) ) ; 

%future  improvement:  add  source  ring  influence 

end 

end 

Vp=Vp+Uvec; 
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A.5  CMV.m 


%Returns  circumferential  mean  tangential , axial ,  and  radial  induced 
%velocities  of  a  propeller  at  any  desired  axial  location 

%Uses  Coney's  version  of  Hough  and  Ordway's  Formulas 

%uHough  and  vHough  are  included  as  a  validation  case.  They  are  only 

%valid  for  no  hub  and  trailing  vortex  system  whose  path  is  determined 
%soley  by  the  incoming  free  stream  with  translation  U  and  rotation 
%(omega).  (Hough  p.319) 


%Variables 

%xC: 

g, 

o 

%rC: 

%M: 

% gamma : 

%rtv: 

%Z: 

%TanBI : 


axl  dist  btwn  prop  (i.e.  vortex  plane)  and  control  point  plane 
positive  if  Ctrl  pt  downstream  of  propeller 
radius  to  calculate  CMV 
#  of  panels  on  lifing  line 
circulation  for  panels  on  lifting  line 

vector  of  radius  of  trailing  vorticies  on  prop  lifting  line 
number  of  blades 

vector  of  tangent  of  advance  angles  of  trailing  vortices 


%Returns :  axlCMV:  axial  CMV  at  xC, rC 

%  radCMV:  radial  CMV  at  xC,  rC 

%  tanCMV:  tangential  CMV  at  xC,  rC 


function  [Vp]  =  CMV(xC, rC, rtv, gamma, TanBI, Z) 


if  abs (xC) <10e-10  &  rC==rtv(end) 

Vp  =  [0;  0;  0] 
return 

end 


%logic  to  skip  routine 
%if  Ctrl  pt  on  lifting  line 


M=length (gamma) ; 

%tangential  velocity  induced  from  a  horseshoe  vortex 

% (bound  and  trailing  vorticity) 

tanCMV=0; 

%  %  %Coney's  implementation.  tanCMV  =  0  if  rtv=rC  or  if  xC=0 
%  %  for  i=l : M 

%  %  S=  (rtv  (i) -rC)  *  (rtv  (i  +  1) -rC)  ; 

%  %  if  S<0  &  xOO 

%  %  tanCMV=tanCMV-Z*gamma (i) / (2*pi*rC) ; 

%  %  end 

%  %  end 

%axial  and  radial  velocity  induced  from  trailing  vorticies 

%sum  effect  from  every  trailing  vortex,  except  for  ends,  each  rtv 
%represents  two  trailing  vortex  with  different  circulation. 

%for  each  horseshoe,  lower  trailer  gamma  is  same  sign  as  bound  gamma 
%and  upper  tailer  is  opposite  sign,  (coney  p.  171) . 
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axlCMV=0 ; radCMV=0 ; 


for  i=l : M+1 

q=l+  (xC''2+  (rC-rtv  (i)  )  ''2)  /  (2*rC*rtv  (i)  )  ; 
s=asin (xC/sqrt (xC^2+ (rC-rtv (i) ) ^2) ) ; 

%amplitude  wrt  elliptical  integrals 
t=sqrt (4*rC*rtv(i) / (xC^2+ ( rC+rtv ( i )  )  ^2)  )  ; 

%t=k;  (modulus  wrt  elliptical  integrals) 

if  rC>rtv(i)  %agrees  with  Coney 

%Hough  has  rc>=rtv(i) 

cl=  xC/ (2*sqrt (rC*rtv (i) ) ) *Q2Mhalf (q) -pi/2*Heuman (s,asin(t) ) ; 

else 

cl=pi+xC/ (2*sqrt (rC*rtv (i) ) ) *Q2Mhalf (q) +pi/2*Heuman (s,asin(t) ) ; 

end 

%  c2  is  not  needed  if  tanCMV  calculated  using  Kelvin's  theorem 
if  rC<rtv(i) 

c2=  xC/ (2*sqrt (rC*rtv (i) ) ) *Q2Mhalf (q) -pi/2*Heuman (s,asin(t) )  ; 

else 

c2=pi+xC/ (2*sqrt (rC*rtv (i) ) ) *Q2Mhalf (q) +pi/2*Heuman (s, asin (t) ) ; 

end 

i~=l  &&  i~=M+l  %logic  for  interior  trailer  vortices 

axlCMV=axlCMV+Z*cl/ (pi*rtv (i)  *TanBI (i) ) * . . . 

(gamma (i) -gamma (i-1) ) ; 

radCMV=radCMV+Z*Q2half (q) / (pi*sqrt (rC*rtv (i) ) *TanBI (i) ) * . . . 

(gamma (i) -gamma (i-1) ) ; 

tanCMV=tanCMV+Z*c2/ (pi*rtv (i)  *  .  .  . 

% (gamma (i) -gamma (i-1) ) ) ; 

else  if  i==l  %logic  for  first  and  last  trailer  vortices 

axlCMV=axlCMV+gamma (i) *Z*cl/ (pi*rtv (i) * . . . 

TanBI (i) ) ; 

radCMV=radCMV+gamma (i) *Z*Q2half (q) / (pi*sqrt (rC*rtv  (i) ) *  .  .  . 
TanBI (i) ) ; 

tanCMV=tanCMV+gamma (i) *Z*c2/ (pi*rtv  (i) ) ; 

else  axlCMV=axlCMV-gamma (i-1) *Z*cl/ (pi*rtv (i) * . . . 

TanBI (i) ) ; 

radCMV=radCMV-gamma (i-1) *Z*Q2half (q) / (pi*sqrt (rC*rtv (i) ) * . . . 
TanBI (i) ) ; 

tanCMV=tanCMV-gamma (i-1) *Z*c2/ (pi*rtv  (i) ) ; 


end 

end 

end 


ax 1 CMV= - ax 1 CMV/ 2 ; 
radCMV=radCMV/2 ; 


%adjust  to  match  PLL  output  for  UA/VS 
%(only  matches  negative  x  values) 
%adjust  to  match  PLL  output  for  UR/VS 


Vp  =  [axlCMV;  radCMV;  tanCMV] ; 
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A.6  ductPlotm 

%Plots  duct 


chordlength 
angle  of  attack 

vortex  ring  radius  (duct  radius  to  meanline) 
chordwise  reference  position  on  duct 
fixed  at  0.5  but  could  be  passed  as  a  variable 
global  propeller  x-coord  of  ductRef 
max  camber  (%  of  chordlength) 
max  thickness  (%  of  chordlength) 

%  %  Notes:  X-axis  positive  in  streamwise  direction  (i.e.  downstream) . 

function []  =  ductPlot (vrRad, c, fo, to, alpha, ductRef ) 

%Read  in  meanline  f (x)  and  thickness  t (x)  distribution  data 
%Read  foil  data  (x,  f/fo,  t/to)  from  text  file 
%Foil_data.txt  contains  parabolic  meanline  (f/fo) 

%and  elliptical  thickness (t/to)  data 


^Variables : 
fe  %  c  [m]  : 

i  %  alpha  [radians] 

i  %  vrRad  [m] : 

fe  %  ductRef: 


xDuct  [m] : 
f  o : 
to : 


[x_over_c, f_over_fo, t_over_to] =textread ( ' foil_data . txt ' ,  ' %f %f %f '  ,  . . . 

' header lines ' , 3 ) ; 

%x  over_c  range  is  -c/2  to  c/2  (this  is  converted  to  x=0  to  x=c  below) 


f=fo*c*f  over  fo; 
t=to*c*t  over  to; 


%camber  distribution 
%thickness  distribution 


x=x  over  c*c  +  c/2; 


%dimensionalizes  x  with  a  range  of  0  to  c 
%range  of  0  to  c  is  needed  for  cosine  spacing 


fpp=spline (x, f ) ; 

%  %  f Ppp=fnder ( fpp) 
tpp=spline (x,  t) ; 


%spline  camber  data 
%splines  slope  of  fpp 
%spline  thickness  data 


%alt  method  not  using  FNDER  (Spline  Toolbox) 
xl=length (x) ; 
theta=zeros (xl, 1) ; 
for  m=l : xl 

if  m==xl 

theta (m) =theta (m-1 ) ; 

else 

theta (m) =atan ( (ppval ( fpp, x (m+1 ) ) -ppval ( fpp, x(m) ) ) / (x(m+l) -x(m) ) ) 

end 

end 


%Generate 
%  %  theta 
x_upper  = 
y_upper  = 
X  lower  = 
y  lower  = 


2-D  flat  cross-section 
=  atan (ppval ( fPpp, x) ) ; 
X  -  t/2 . *sin (theta) ; 
f  +  t/2 . *cos (theta) ; 

X  +  t/2 . *sin (theta) ; 
f  -  t/2 . *cos (theta) ; 
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%  %  %  %Plot  2-D  section  with  0  degrees  angle  of  attack 
%  %  %  plot (x_lower , y_lower ) 

%  %  %  hold  on 

%  %  %  plot (x_upper , y_upper ) 

%  %  %  xlabel (' X-axis '); ylabel (' Y-axis ') ; 

%  %  %  title ('Duct  section  with  0  degrees  angle  of  attack'); 
%  %  %  axis  equal 
%  %  %  figure 


%Reposition  section 

varl=ductRef *c;  %offset  (x-dir) ,  ductRef  at  x=0 

var2=ppval (fpp, ductRef *c) ;  %offset  (y-dir)  for  camber,  ductRef  at  y=0 

var3=0 . 5*ppval (tpp, ductRef *c) ;  %offset  (y-dir)  for  thick,  ductRef  at  y=0 

%ductRef  is  point  where  blade  and  duct  meet 

X  upper  =  X  upper  -  varl; 

y  upper  =  y  upper  -  var2  +  var3; 

X  lower  =  X  lower  -  varl; 

y  lower  =  y  lower  -  var2  +  var3; 


%Rotate  for  angle  of  attack  and  place  section  at  correct  radius 
X  upper  rot  =  x  upper*cos ( -alpha)  -  y  upper*sin (-alpha) ; 

y  upper  rot  =  x  upper*sin (-alpha)  +  y  upper*cos ( -alpha)  +  vrRad; 

X  lower  rot  =  x  lower*cos ( -alpha)  -  y  lower*sin (-alpha) ; 

y  lower  rot  =  x  lower*sin (-alpha)  +  y  lower*cos ( -alpha)  +  vrRad; 

%Plot  duct  section  rotated 
plot (x_lower_rot, y_lower_rot) 
hold  on 

plot (x_upper_rot, y_upper_rot) 
xlabel ( ' X-axis ' ) ; ylabel ( ' Y-axis ' ) ; 

title (['Duct  section  (repositioned)  with  ' , num2str (alpha*180/pi) , . . . 

'  degrees  angle  of  attack']); 

axis  equal 
figure 


5-9-9- 
o  o  o 

9-9-9- 
o  o  o 

9-9-9- 
o  o  o 

9-9-9- 
o  o  o 

9-9-9- 
o  o  o 

9-9-9- 
o  o  o 

9-9-9- 
o  o  o 

9-9-9- 
o  o  o 

9-9-9- 
o  o  o 


%Build  all  sections  (upper  and  lower  surfaces)  for  complete  3-D  duct 
z= zeros (length (x) , 1 ) ; 

[thetaU, phiU, RU] =cart2sph (y_upper_rot, x_upper_rot, z) ; 

[thetaL, phiL, RL] =cart2sph (y  lower  rot,x  lower  rot,z); 


nds=50;  %#  of  duct  sections  for  plotting 

for  n=l : nds 

%  phi=0+pi : 1 . 9*pi/ (nds-1 ) : 2 . 1 *pi+pi ;  %360  deg  coverage  for  duct 

phi=0 : 2 *pi/ (nds-1 ) : 2 . l*pi ;  %360  deg  coverage  for  duct 

[x_u_3D (n, : ) , y_u_3D (n, : ) , z_u_3D (n, : ) ] =sph2cart (phi (n) , thetaU, RU) ; 
[x_l_3D (n, : ) , y_l_3D (n, : ) , z_l_3D (n, : ) ] =sph2cart (phi (n) , thetaL, RL) ; 


%  %  %  %Plot  duct  sections  individually 

%  %  %  plot3 ( z_u_3D (n, : ) , x_u_3D (n, : ) , y_u_3D (n, : ) ) 

%  %  %  hold  on 

%  %  %  plot3 ( z_l_3D (n, : ) , x_l_3D (n, : ) , y_l_3D (n, : ) ) 

end 

%  %  %  xlabel (' X-axis '); ylabel (' Y-axis '); zlabel (' Z-axis '  ) 

%  %  %  title (['Duct  with  ' , num2str (alpha*180/pi) , '  degrees  angle  of  attack']) 
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o\o  o\o 


%  %  %  axis  equal 
%  %  %  figure 

%Plot  duct  as  3-D  surface 
surfl(z  u  3D,x  u  3D,y  u  3D) 
hold  on 

surfl(z  1  3D,x  1  3D,y  1  3D) 

%  %  %  xlabel (' X-axis '); ylabel (' Y-axis '); zlabel (' Z-axis ' ) 

%  %  %  title (['Duct  with  ' , num2str (alpha*180/pi) , '  degrees  angle  of  attack']) 
%  %  %  axis  equal 

shading  interp 
colormap (copper) 
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Appendix  B.  Mathematical  Functions  MATLAB®  Code 


B.l  Qlhalf.m 

%Q2half:  Legendre  fuction  of  the  second  kind  and  positive  half  order 
%Ref:  Handbook  of  Math  Functions,  Abramowitz  and  Stegun,  1972 
%section  8.13.7,  p.337 

%uses  modulus  k  for  elliptic  integrals  (m=k^2) 

%ellipke  uses  parameter  m. 

function  [Q2]  =  Q2half (q) 

k=sqrt (2/ (q+1) ) ; 

[K, E] =ellipke (k^2) ; 

Q2=q*k*K-sqrt (2* (q+1) ) *E; 

%Validated  with  the  National  Bureau  of  Standards  Tables  of 
%Associated  Legendre  Functions 

% (Columbia  University  Press,  New  York,  1945),  p.266. 

%From  the  tables:  Q2half (1 . 5) =. 393175,  Q2half (2 . 7) =. 134035, 

%Q2half (6) =.  0382887,  Q2half (8 . 4) =. 0229646,  Q2half (10) =. 0176449 


B.l  QlMhalf.m 

%Q2Mhalf:  Legendre  fuction  of  the  second  kind  and  minus  half  order 

%Ref:  Handbook  of  Math  Functions,  Abramowitz  and  Stegun,  1972 
%section  8.13.3,  p.337 

%uses  modulus  k  for  elliptic  integrals  (m=k^2) 

%ellipke  uses  parameter  m. 

function  [Q2M]  =  Q2Mhalf (q) 

k=sqrt (2/(q+l)); 

[K, E] =ellipke (k^2) ; 

Q2M=k*K; 


%checked  with  ref  p.340  example 

%Validated  with  the  National  Bureau  of  Standards  Tables  of 
%Associated  Legendre  Functions 

% (Columbia  University  Press,  New  York,  1945),  p.264. 

%%From  the  tables:  Q2Mhalf (1 . 5) =2 . 01891,  Q2Mhalf (2 . 7) =1 . 38958, 
%Q2Mhalf (6) =0.911696,  Q2Mhalf (8 . 4) =0 . 768523,  Q2Mhalf (10) =0 . 703806 
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o\°  o\°  o\°  o\°  o\° 


B.3  Heuman.m 


Heuman:  Heuman ' s  Lambda  faction 

Ref:  Handbook  of  Math  Functions,  Abramowitz  and  Stegun,  1972 
%section  17.4.39,  p.595 

Ref:  Handbook  of  Elliptic  Integrals  for  Engineers  and  Physicists,  Byrd  and 
%Friedman,  1954,  p37. 

phi:  amplitude  (radians),  (CMV  sends  's'  as  phi) 

alpha:  modular  angle  (radians),  (CMV  sends  't'  alpha) 

function  [H]  =  Heuman (phi , alpha) 

[K,  E]  =ellipke  (sin  (alpha)  ''2)  ; 

F=mfun ( ' EllipticF ' , sin (phi ) , sin (pi/ 2 -alpha) ) ; 

%lncomplete  elliptic  integral,  1st  kind 
EE=mfun ( ' EllipticE ' , sin (phi ) , sin (pi/ 2 -alpha) ) ; 

%Incomplete  elliptic  Integral,  2nd  kind 


H=2/pi* (K*EE- (K-E) *F) ; 
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Appendix  C.  Variational  Optimization  Routine  MATLAB®  Code 


C.l  Coney. m 


g, _ 

o - 

%  ==========================================================  Coney  Function 

0, 

o 

%  The  Coney  function  determines  the  "optimum"  circulation  distribution 
%  that  satisfies  the  input  operating  conditions,  using  a  variational 
%  optimization  algorithm,  as  described  on  Coney,  page  25.  The  Coney 
%  function  returns  performance  specs,  such  as  thrust  coefficient  and 
%  efficiency,  as  well  as  the  circulation  distribution,  ect. 

g, 

o 

%  Reference:  Coney,  William,  "A  Method  for  the  Design  of  a  Class  of  Optimum 
%  Marine  Propulsors",  Ph . D .  thesis,  MIT,  1989. 

g, 

o 

%  Includes  Coney  and  Align  wak  functions 

%  Authors:  Brenden  Epps  (variational  optimization  and  wake  alignment) 

%  Mitch  Stubblefield  (duct  integration) 

g, _ 

o 


function  [CT,  CQ,  CP,  KT, KQ, VMIV, EFFY, RC, G, VAC, VTC, UASTAR, UTSTAR, TANBC,  .  .  . 

TANBIC, CoD, CD, TAU, Xring, dVort, UADUCT, dCirc, UAdVS,  URdVS]  . .  . 

=  Coney (Rhub, R, Z,Mp, ITER, Rhv, HUF, TUF, SCF, Js , CTDES , Hub_Flag, . . . 
Duct_Flag, TAU, rDuct_oR, CDd, XR, XCoD, XCD, XVA, XVT, rho, Vs) ; 

clc 


rDuct=rDuct  oR*R; 
UADUCT=zeros (l,Mp)  ; 
CTD=0; 

dVort=zeros ( 1 , Mp+2 ) ; 


Initialize  variables  needed  in  functions 
%duct  radius 

%Induction  of  duct  on  lifting  line  Ctrl  pts 
%CT  for  duct 

%circulation  distribution  of  each  vortex  ring 


Xring=zeros (l,Mp+2) ; UAdVS=Xring; URdVS=Xring; 


% -  Compute  the  Volumetric  Mean  Inflow  Velocity,  eqn  163,  p.l38 

Rhub  oR  =  Rhub/R;  %  [  ] ,  hub  radius  /  propeller  radius 

RoR  =1;  %  [  ],  propeller  radius  /  propeller  radius 


VMIV  =  2*trapz  (XR,  XR.  *XVA)  /  (RoR''2-Rhub_oR''2)  ;  %  [  ],  VMIV/ship  velocity 


%  -  Compute 

RV=zeros ( 1 , Mp+1 ) ; RC=zeros ( 1 , Mp) ; 
if  Duct  Flag==0  &  Hub  Flag==0 

DRR~=  (RoR-Rhub_oR) / (Mp+ .5); 

RV ( Mp+ 1 ) =RoR- .25*  DRR; 

RV(1) =Rhub_oR+.25*DRR; 
elseif  Duct  Flag==l  &  Hub  Flag==0 
DRR  =  (RoR-Rhub_oR)/(Mp+.25); 
RV(Mp+l) =RoR; 

RV(l)  =Rhub_oR+.25*DRR; 
elseif  Duct  Flag==l  &  Hub  Flag==l 
DRR  =  (RoR-Rhub_oR) / (Mp) ; 
RV(Mp+l) =RoR; 

RV(l)=Rhub  oR; 


evenly-spaced  vortex  &  control  pt .  radii 

%  initialize  RC  and  RV 
%  no  duct  image  or  hub  image 
%  panel  size 
%  25%  tip  inset 
%  25%  hub  inset 
%  duct  image  but  no  hub  image 
%  panel  size 
%  no  tip  inset 
%  25%  hub  inset 
%  duct  image  and  hub  image 
%  panel  size 
%  no  tip  inset 
%  no  hub  inset 
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elseif  Duct  Flag==0  &  Hub  Flag==l 
DRR  =  (RoR-Rhub_oR)/(Mp+.25); 
RV ( Mp+ 1 ) =RoR- .25*  DRR; 

RV(1) =Rhub_oR; 

end 

RC(1)=RV(1)+.5*DRR; 
for  m=2 : Mp 

RV(m) =RV(m-l) +DRR; 

RC  (m)  =RC  (m-1)  +DRR; 

end 


%  no  duct  image  but  hub  image 
%  panel  size 
%  25%  inset  for  tip 
%  no  hub  inset 


%  Ctrl  pt  at  mid-panel 


DR  =  diff (RV) ; 


%  difference  in  vortex  radii  /  propeller  radius 


%  -  Interpolate 

VAV  =  pchip (XR, XVA, RV) ; 

VTV  =  pchip (XR, XVT, RV) ; 

VAC  =  pchip  (XR,  XVA,  RC)  ; 

VTC  =  pchip (XR, XVT, RC) ; 

CD  =  pchip (XR, XCD, RC) ; 

CoD  =  pchip (XR,XCoD,RC); 


Va,  Vt,  Cd,  and  c/D  at  vortices  &  control  points 
%  axial  inflow  vel.  /  ship  vel.  at  vort  pts 

%  tangential  inflow  vel.  /  ship  vel.  at  vort  pts 

%  axial  inflow  vel.  /  ship  vel.  at  Ctrl  pts 

%  tangential  inflow  vel.  /  ship  vel.  at  Ctrl  pts 

%  section  drag  coefficient  at  Ctrl  pts 

%  section  chord  /  propeller  diameter  at  Ctrl  pts 


TANBC  =  VAC./ (pi.*RC./Js  +  VTC) ; 


%  tan (Beta)  at  control  pts. 


%Allocate  CTDES  between  propeller  and  duct 

CTPDES  =  CTDES*TAU;  %CT  desired  for  the  propeller 

CTDDES  =  CTPDES/TAU-CTPDES;  %CT  desired  for  the  duct 

VD  =0;  %viscous  drag 

%Initial  guess  for  dCirc  (circulation  on  duct) 
dCirc  =  0.5* (1-TAU) ; 

if  TAU==1  &  CDd~=0  %provides  a  small  duct  circulation  to  offset  drag 

dCirc  =  .001; 

end 


%  -  Compute  vortex  ring  influence  functions  from  duct 

if  Duct  Flag  ==  1 

[vRingLoc, dVort, UADUCT]  =  ductVort (R, rDuct, Mp, RC) ; 
Xring=vRingLoc (1, ; ) ; 

end 


%  =====================  Determine  optimum  circulation  distribution  function 

g, 

o 

%  See  William  Coney  Ph.D.  thesis  "A  Method  For  the  Design  of  a  Class  of 
%  Optimum  Marine  Propulsors",  MIT,  1989,  page  25  for  a  discussion  of 
%  this  variational  optimization  algorithm. 


%  -  Initialize 

UASTAR(l:Mp)  =  0; 

UTSTAR(l:Mp)  =  0; 

LM  last  =  -1; 

G  last  =  0; 


induced  velocities  &  Lagrange  multiplier  (Coney  p.27) 

%  UASTAR  /  ship  speed 
%  UTSTAR  /  ship  speed 

%  last  value  of  the  Lagrange  Multiplier, 
%  last  value  of  the  circulations,  G 


LM 


A  =  zeros (Mp+1); 

B  =  zeros (Mp+1 , 1 ) ; 


%  A  matrix  for  linear  system  of  equations 
%  B  matrix  for  linear  system  of  equations 
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%Estimate  the  axial  induced  velocity  with  actuator  disc  approx 
for  m=l ; Mp 

UASTAR(m) =0 . 5* (sqrt (1+CTPDES) -1)  +  dCirc*UADUCT (m) ; 

end 

%Initial  TANBIC  and  TANBIV  estimates 

[TANBIC, TANBIV]  =  f ind_tan_BetaI (VAC, VTC, UASTAR, UTSTAR, RC, RV, Js) ; 

%  -  Iterate  to  solve  simultaneous  equations  for  G,  LM,  &  Betal .  Fix 

%  Betal,  and  solve  simultaneous  equations  for  G  and  LM. 

B  iter  =  1;  %  iteration  in  the  Betal  loop 

B  res  =1;  %  residual  Betal  between  interations 

TANBIC  last  =  TANBIC;  %  the  last  value  of  TANBIC 


while  B_iter  <  ITER  &  B_res  >  Ie-5  %  (WHILE  LOOP  Bl) 

%  -  Compute  the  vortex  Horseshoe  Influence  Functions,  p.l79 

[UAHIF, UTHIF]  =  Horseshoe (Mp, Z , TANBIV, RC , RV, SCF, Hub_Flag, Rhub_oR, ... 

Duct_Flag, rDuct_oR) ; 

%  -  Iterate  to  solve  simultaneous  equations  for  G  and  LM  for  the 

%  current  values  of  Betal.  (Coney  eqns .  2.32  &  2.33,  p.27) 

G  iter  =1;  %  iteration  in  the  G  loop 

G  res  =1;  %  residual  G  between  interations 

LM  res  =1;  %  residual  LM  between  interations 

while  G_iter  <  ITER  &  (G_res  >  le-5  |  LM_res  >  le-5)  %  (WHILE  LOOP  GI) 

%  -  Solve  simultaneous  equations  for  G  and  LM 

for  i  =  l:Mp  %  for  each  vortex  panel,  i 

for  m  =  l;Mp  %  for  each  vortex  panel,  m 

A(i,m)  =  UAHIF (i,m) *RC (m) *DR(m)  +  ...  %  A 

UAHIF (m, i) *RC (i) *DR(i)  +  ... 

LM_last*UTHIF (i,m) *DR(m)  +  ... 

LM_last*UTHIF (m, i) *DR(i)  ; 

end 

A(i,Mp+l)  =  (VTC(i)  +  pi*RC(i)/Js)  *DR(i);  %  C 

A(Mp+l,i)  =  (VTC(i)  +  pi*RC(i)/Js  +  UTSTAR ( i )) *DR ( i ) ;  %  D 

B(i)  =  - (VAC (i) +dCirc*UADUCT (i) ) *RC (i) *DR(i) ;  %  B 

end 

B(Mp+l)  =  (CTPDES+VD) / (4*Z) ;  %  D 

GLM  =  linsolve (A, B) ;  %  Solve  Mp+1  by  Mp+1  system  of  equations 

G  =  GLM(l:Mp);  %  G  is  the  first  Mp  entries 

LM  =  GLM(Mp+l);  %  LM  is  the  last  entry 

%  -  Compute  induced  velocities  at  control  points  eqn  254,  p.l79 
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[UASTAR, UTSTAR]  =  Induced_Velocity (Mp, G, UAHIF, UTHIF, UADUCT, dCirc) ; 


%  -  Calculate  duct  thrust  (including  propeller  influence) 

%  and  scale  duct  circulation  (dCirc) 

if  Duct  Flag  ==  1 

[CTD, dCirc, UAdVS, URdVS]  =  ductThrust (vRingLoc, dVort, dCirc, .. . 

rDuct, CDd, XVA(end) , rho, RV, C , TANBIV, Z, R, CTDDES) ; 


end 


% -  Prepare  for  the  next  iteration 

G  iter  =  G  iter  +1  %  iteration  in  the  G  loop 

G_res  =  abs (G  -  G_last) ;  %  residual  G  between  interations 

LM  res  =  abs (LM  -  LM  last) ;  %  residual  LM  between  interations 

G_last  =  G;  %  the  last  value  of  G 

LM  last  =  LM;  %  the  last  value  of  LM 

if  G_iter  >  ITER 
warning ( ' on ' ) , 

warning (' WARNING :  While  loop  G1  did  NOT  converge.'), 
warning ( ' of f ' ) , 

end 

end  %  (END  WHILE  LOOP  GI) 

% - Align  the  wake  to  the  new  circulation  distribution 

[UAHIF, UTHIF, UASTAR, UTSTAR, TANBIC, TANBIV]  =  ... 

Align_wake (TANBIC, TANBIV, ITER,Mp, Z, RC, RV, SCF, Hub_Flag, Rhub_oR, . . . 

G, VAC, VTC, Js, Duct_Flag, rDuct_oR, UADUCT, dCirc) ; 

%  -  Prepare  for  the  next  iteration 

B  iter  =  B  iter  +1  %  iteration  in  the  Betal  loop 

B  res  =  abs (TANBIC  -  TANBIC  last)  ;  %  residual  Betal  between  interations 
TANBIC_last  =  TANBIC;  %  the  last  value  of  TANBIC 

if  B_iter  >  ITER 
warning ( ' on ' ) , 

warning (' WARNING :  While  loop  B1  did  NOT  converge.'), 
warning ( ' of f ' ) , 

end 

[CT, CQ, CP, KT, KQ, EFFY, TAU, VD]  =  Forces (CD, RV, VAC, TANBC, UASTAR, UTSTAR, ... 

CoD, G,Mp, RC, Hub_Flag, Rhv, Z, Js, VMIV, CTDDES) ; 

end  %  (END  WHILE  LOOP  Bl) 

%  -  Compute  thrust  &  torque  coefficients  and  efficiency 

[CT, CQ, CP, KT, KQ, EFFY, TAU, VD]  =  Forces (CD, RV, VAC, TANBC, UASTAR, UTSTAR, . . . 

CoD, G,Mp, RC, Hub_Flag, Rhv, Z, Js, VMIV, CTDDES) ; 

%  %  -  If  required,  unload  the  hub  and  tip,  then  rescale  the  circulation 

%  %  distribution  to  get  the  desired  value  of  the  thrust  coefficient. 

%  if  Hub_Flag  &  (HUF  >  0  |  TUF  >0)  %  (IF  STATEMENT  Ul) 

g, 

o 

%  %  - Unload  hub  and  tip  as  specified  by  HUF  and  TUF 

%  RU  =  (RC  -  Rhub  oR)./(RoR  -  Rhub  oR) ; 
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g, 

o 

g, 

o 

q, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 


nH  =  4; 
nT  =  3; 

GH  =  HUF*G(1)  *sqrt  (1-RU. ''2)  .  *  (1-RU.  ^'2)  . (2*nH-2)  ; 

GT  =  TUF*G  (Mp)  *sqrt  (1-RU.  ■^2)  .  *  (  RU  .  ^^2  )  .  ^  ( 2  *nT-2  )  ; 

G  =  G  -  GH '  -  GT '  ; 

%  =============  END  determine  optimum  circulation  distribution  function 

% - Align  the  wake  to  the  new  circulation  distribution 

[UAHIF, UTHIF, UASTAR, UTSTAR, TANBIC, TANBIV]  =  ... 


Align_wake (TANBIC, TANBIV, ITER,Mp, Z, RC, RV, SCF, Hub_Flag, Rhub 


oR,  G, VAC, VTC,  Js,  .  . 


g, 

o 

g, 

o 

g, 

o 

g, 

o 

q, 

o 

g, 

o 

q, 

o 

q, 

o 

g, 

o 

q, 

o 

g, 

o 

q, 

o 

q, 

o 

g, 

o 

q, 

o 

g, 

o 

q, 

o 

q, 

o 

g, 

o 

q, 

o 

g, 

o 

q, 

o 

g, 

o 

g, 

o 

q, 

o 

g, 

o 

q, 

o 

g, 

o 

q, 

o 

q, 

o 

g, 

o 


Duct  Flag, rDuct  oR, UADUCT, dCirc) 


CT_iter 

CT_res 

CT_last2 

CT_lastl 

GMF_last2 

GMF  lastl 


Iterate  to  scale  G  to  get  desired  value  of  thrust  coefficient 
=  1 ;  %  iteration  in  the  CT  loop 

=  1;  %  residual  CT 

=  0;  %  the  CT  prior  to  the  last  CT 

=  0;  %  the  last  value  of  CT 

=  0;  %  the  GMF  prior  to  the  last  GMF 

=  0;  %  the  last  value  of  GMF 


while  CT  iter  <  ITER  &  CT  res  >  Ie-5 


(WHILE  LOOP  CTI) 


if  CT_iter  ==  1 
GMF  =  1; 


elseif  CT_iter  ==  2 

GMF  =  1+ (CTPDES-CT) / (5*CTPDES) ; 


elseif  CT  iter  >  2 

GMF  =~GMF_lastl  +  (GMF_lastl-GMF_last2 ) * (CTPDES-CT_lastl ) / . . . 

(  CT_lastl-  CT_last2); 

end 


G  =  GMF.*G;  %  GMF  =  G  Multiplication  Factor 

%  -  Compute  induced  velocities  at  control  points,  eqn  254,  p.l79 

[UASTAR, UTSTAR]  =  Induced_Velocity (Mp, G, UAHIF, UTHIF, UADUCT, dCirc) ; 

%  -  Compute  thrust  &  torque  coefficients  and  efficiency 

[CT, CQ, CP, KT, KQ, EFFY, TAU, VD]  =  ... 


Forces (CD, RV, VAC, TANBC, UASTAR, UTSTAR, CoD, G,Mp, RC, Hub  Flag, Rhv, Z, Js, VMIV, CTD) ; 


g, 

o 

g, 

o 

q, 

o 

g, 

o 

q, 

o 

g, 

o 

g, 

o 

q, 

o 

g, 

o 

q, 

o 


q,  _  _ 

o 

CT_iter 

CT_res 

CT_last2 

CT_lastl 

GMF_last2 

GMF  lastl 


=  CT_iter  +  1 ; 

=  abs (CT  -  CTPDES 
=  CT_lastl; 

=  CT; 

=  GMF_lastl; 

=  GMF; 


-  Prepare  for  the  next  iteration 
%  iteration  in  the  CT  loop 
%  residual  CT 

%  the  CT  prior  to  the  last  CT 
%  the  last  value  of  CT 
%  the  GMF  prior  to  the  last  GMF 
%  the  last  value  of  GMF 


end 


%  (END  WHILE  LOOP  CTI) 
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o\°  o\°  o\°  o\° 


end 


(END  IF  STATEMENT  Ul) 


END  Coney  Function 


q, 

o 

g, 

o 

0, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

q, 

o 

g, 

o 

q, 

o 


=====================================================  Align  wake  Function 

This  function  aligns  the  wake  to  the  given  circulation  distribution  by 
iteratively  computing: 

UAHIF,UTHIF  =  the  horseshoe  influence  functions 
UASTAR, UTSTAR  =  the  induced  velocities 
TANBIC, TANBIV  =  the  velocity  angles 
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o\o  o\o 


C.2  AUgnwake.m 


function  [UAHIF, UTHIF, UASTAR, UTSTAR, TANBIC, TANBIV]  =  ... 

Align_wak;e  (TANBIC,  TANBIV,  ITER,Mp,  Z,  RC,  RV,  SCF,  Hub_Flag,  Rhub_oR,  .  .  . 
G, VAC, VTC, Js, Duct_Flag, rDuct_oR, UADUCT, dCirc) 

%  -  Iterate  to  ALIGN  WAKE  to  the  new  circulation  distribution 

W  iter  =1;  %  iteration  in  the  wake  alignment  loop 

W  res  =1;  %  residual  Betal  between  interations 

TANBIW_last  =  TANBIC;  %  the  last  value  of  TANBIC 

while  W_iter  <  ITER  &  W_res  >  Ie-5  %  (WHILE  LOOP  WAl) 

%  -  Compute  the  vortex  Horseshoe  Influence  Functions,  p.l79 

[UAHIF, UTHIF]  =  Horseshoe (Mp, Z , TANBIV, RC , RV, SCF, Hub_Flag, .. . 

Rhub  oR, Duct  Flag, rDuct  oR) ; 

%  -  Compute  induced  velocities  at  control  points,  eqn  254,  p.l79 

[UASTAR, UTSTAR]  =  Induced_Velocity (Mp, G, UAHIF, UTHIF, UADUCT, dCirc) ; 

%  -  Compute  tan (Betal)  for  the  new  induced  velocities 

[TANBIC, TANBIV]  =  f ind_tan_BetaI (VAC, VTC, UASTAR, UTSTAR, RC, RV, Js) ; 

% -  Prepare  for  the  next  iteration 

W  iter  =  W  iter  +1  %  iteration  in  the  Betal  loop 

W  res  =  abs (TANBIC  -  TANBIW  last) ;  %  residual  Betal 
TANBIW_last  =  TANBIC;  ~  %  the  last  value  of  TANBIC 

if  W_iter  >  ITER 
warning ( ' on ' ) , 
warning ( ' WARNING : 
warning ( ' of f ' ) , 

end 

end 


END  Align  wake  Function 


While  loop  WAl  did  NOT  converge.'), 

%  (END  WHILE  LOOP  WAl) 
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Appendix  D.  Test  Case  Setup  Data 


PLL 

OpenProp  v2 

Advance  coefficient,  Js 

0.60 

Thrust  Coefficient,  Cj 

1.20 

Blade  number 

5 

Speed  (propeller) 

150  RPM 

Drag  Coefficient 

0.008 

Vortex  panels  (Mp) 

10 

Diameter  (prop) 

10  ft 

3.048  m 

Diameter  (duct) 

10  ft 

3.048  m 

Speed  (ship) 

15  ft/s 

4.572  m/s 

Thrust 

21206  Ibf 

94328  N 

Figure  D-1:  Test  case  parameters 


OpenProp 


Options 


1 

5 

Number  of  Blades 

l~  Hub  Image  Flag  (Check  for  YES) 

n- 

Thrust  Ratio 

r 

Propeller  Speed  (RPM) 

17  Ducted  propeller  (Check  for  YES) 

1  1 

Duct  Diameter/Prop  Diameter 

I  0.008| 

)uct  Section  Drag  Coefficient 

1 

3.048 

Propeller  Diameter  (m) 

Meanline  Type: 

Thickness  Form: 

|naca  a=0.8 

“3 

|HACA65A010 

J 

1 

94328 

Required  Thrust  (N) 

1  r/R 

1  1 

Cd 

1 

VaA/s  I 

VWs 

fO/c  I 

tO/c  1 

Skew 

Xs/D 

r 

4.572 

Ship  Velocity  (m/s) 

1  0.2 

1  0.16  1 

0.008 

T 

1  1 

0 

0.0174  1 

0.2056  1 

0 

0 

r 

0.6096 

Hub  Diameter  (m) 

1  0.3 

1  0.1818  1 

0.008 

1 

1  1 

0 

0.0195  1 

0.1551  1 

0 

0 

r 

10 

Number  of  Vortex  Panels  over  the  Radius 

1  0.4 

1  0.2024  1 

0.008 

T 

1  1 

0 

0.0192  1 

0.1181  1 

0 

0 

r 

10 

Max.  Iterations  in  Wake  Alignment 

1  0.5 

1  0.2196  1 

0.008 

T 

1  1 

0 

0.0175  1 

0.0902  1 

0 

0 

r 

1 

Hub  Vortex  Radius/Hub  Radius 

1  0.6 

1  0.2305  1 

0.008 

1 

1  1 

0 

0.0158  1 

0.0694  1 

0 

0 

r 

0 

Hub  Unloading  Factor:  0=Optimum 

1  0.7 

1  0.2311  1 

0.008 

1 

1  1 

0 

0.0143  1 

0.0541  1 

0 

0 

r 

0 

Tip  Unloading  Factotr:  1=Reduced  Loading 

1  0.8 

1  0.2173  1 

0.008 

T 

1  1 

0 

0.0133  1 

0.0419  1 

0 

0 

r 

1 

Swirl  Cancellation  Factor:  1=No  Cancellation 

1  0.9 

1  0.1806  1 

0.008 

T 

1  1 

0 

0.0125  1 

0.0332  1 

0 

0 

r 

1031 

Water  Density  (l(g/m*3) 

1  0.95 

1  0.1387  1 

0.008 

T 

1  1 

0 

0.0115  1 

0.0324  1 

0 

0 

1  ^ 

1  0.001  1 

0.008 

T 

1  [ 

0 

0  1 

0  1 

0 

0 

r 

3.048 

Shaft  Centerline  Depth  (m) 

Filename  Prefix 

r 

0.3 

Inflow  Variation  (m/s) 

1 

OpenProp 

r 

1.54 

Ideal  Angle  of  Attack  (degrees) 

1 

Number  of  Points  over  the  Chord 

Run  OpenProp 

Figure  D-2:  OpeuProp  v2  iuput  GUI  for  test  cases 
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PLL  CURRENT  SETTINGS  (inviscid  runs) 


1... 

..circulation  optimum  enabled. ...T 

2... 

.chord  optimization  enabled. 

. F 

3... 

..wake  alignment  disabled 

. F 

4... 

.a  =  0.8  meanline  for  duct . 

. T 

5... 

..number  of  panels . 

. 10 

6... 

.hub  vortex  radius . 

.  1.00 

7... 

..tip  thickness/diameter... 

0.0040 

8... 

.Lagrange  multiplier . 

-1.00 

9... 

..max.  lift  coefficient . 

0.6000 

10. 

...max.  thickness/chord . 

.  0.20 

11.. 

...minimum  root  chord . 

0.1600 

12. 

....enable  computed  drag  coeff . T 

13.. 

...drag  coeff.  multiplier . 

0.0000 

14. 

....duct  tip  gap  factor . 

.  1.00 

PLL  CURRENT  SETTINGS  (viscid  runs) 


1... 

..circulation  optimum  enabled. ...T 

2... 

.chord  optimization  enabled. 

. F 

3... 

..wake  alignment  disabled 

. F 

4... 

.a  =  0.8  meanline  for  duct . 

. T 

5... 

..number  of  panels . 

. 10 

6... 

.hub  vortex  radius . 

.  1.00 

7... 

..tip  thickness/diameter... 

0.0040 

8... 

.Lagrange  multiplier . 

-1.00 

9... 

..max.  lift  coefficient . 

0.6000 

10. 

...max.  thickness/chord . 

,  0.20 

11.. 

...minimum  root  chord . 

0.1600 

12. 

....enable  computed  drag  coeff . F 

13.. 

...drag  coeff.  multiplier . 

0.0080 

14. 

....duct  tip  gap  factor . 

.  1.00 

Figure  D-3:  PLL  current  settings  for  inviscid  and  viscid  test  cases 


PROPELLER  LIETING  LINE  RUN:  eta_caEe  00/00/94 

OVERALL  INPUT  FILE 


15.0000  .  ship  speed  (ft/sec) 

2.0000  .  Fluid  density 

10.0000  .  Shaft  centerline  depth  (ftO 

1  Number  of  components 

N  No  image  hub  to  be  used 

Y  Image  duct  to  be  used 

0.5000  .  (Duct  chord  1 ength)/ (Component  #1  diameter) 

0.0000  .  Drag  coefficient  for  the  duct 

0.0000  .  (Duct  thi ckness)/(Component  #1  diameter) 

10.0000  .  Duct  diameter  (ft) 

0.0000  .  Axial  location  of  duct  mid-chord  (ft) 

5  Number  of  blades  on  component  1 

10.0000  .  Diameter  of  component  1  (ft) 

a. bid  File  containing  blade  inputs  for  comp.  1 

10.0000  .  Diameter  of  wake  for  component  1 

a.wak  File  containing  wake  inputs  for  component  1 


Figure  D-4:  PLL  overall  input  file  for  test  cases 
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PROPELLER  LIETING  LINE  RUN:  eta_case  00/00/94  00:0 

Js=0 . 6  ,  Ct=l .  2  ,  Invi  send 

PLL  VERSION  4.2  SUMMARY  OUTPUT 


SHIP  SPEED  Oft/sec):  15.00 
SHAET  DEPTH  (ft):  10.00 
TOTAL  THRUST  OlbsO:  21203.93 
EPPICIENCY:  0.764 
IMAGE  HUE  USED:  E 


WAKE  ALIGN.  DISABLED: 
TIP  GAP  PACTOR: 

DUCT  DIAMETER 


PLUID  DENSITY:  2.0000 
#  PROP.  COMPONENTS:  1 
HORSEPOWER:  756.539 
VOL.MN.VEL. (l-W):  1.000 
IMAGE  DUCT  USED:  T 
CIRCULATION  OPTIMIZED:  T 
CHORDLENGTHS  OPTIMIZED:  E 
DUCT  AXIAL  LOC.  fft/ :  0.000 


CftO: 


E 

1.0000 

10.000 


DUCT  OUTPUTS 


DUCT  AREA  RATIO: 

2.000 

DUCT 

VOLUME  Ccu  ft): 

0.00 

DUCT  THRUST  (lbs) : 

0.64 

DUCT 

DIAMETER  (ft): 

10.0000 

TAU  CTp/TO: 

1.0000 

DUCT 

KT: 

0.0000 

CHD/PROP  DIAM: 

0.5000 

DRAG 

COEFFICIENT: 

0.0000 

THICKNESS/P  DIAM: 

0.0000 

DUCT 

CIRC. : 

0.0000 

CT: 

0.0000 

CL: 

0.0001 

AXIAL  INFLOW  VEL : 

1.0000 

RADIAL  INFLOW  VEL: 

0.0000 

VISCOUS  DRAG  (lbs): 

0.00 

INV. 

THRUST  (lbs): 

0.64 

SUMMARY  OUTPUT  POR  COMPONENT  NUMBER  1 


AXIAL  LOCATION (ft) : 

0.00 

NUMBER  OP  BLADES: 

5 

NUMBER  OF  PANELS: 

BLADE  PILE:  a. bid 

WAKE  FILE:  a.wak 

10 

TORQUE  (ft-lbs): 

26489.53 

3  SHIP: 

0.600 

BLADE  VOLUME  Cft**3): 

23.1175 

EXPANDED  AREA  RATIO: 

0.  508 

DIAMETER  (ft):  10.00 
HUE  DIAMETER  (ft):  2.00 
REVS.  PER  MINUTE:  150.00 
WAKE  DIAMETER  (ft):  10.00 
THRUST  tlbsO:  21203.29 
HORSEPOWER:  756.539 
VOLUMETRIC  J:  0.600 
MOM  INERTIA  Cft**5):  157.17 


CT:  1.1999  CQ:  0.1499  CP:  1. 5697  KT:  0.1696  KQ:  0.0212 


Figure  D-5:  Sample  PLL  output  summary  for  test  case  ruu 
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